text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMDynamics.jl/blob/master/LICENSE

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess
using Logging
Logging.configure(level=INFO)
add_elements! = JuliaFEM.add_elements!
using FEMDynamics

# Model construction starts

datadir = Pkg.dir("FEMDynamics", "examples", "ball")
meshfile = joinpath(datadir, "ball.med")
mesh = aster_read_mesh(meshfile)
mesh.element_sets[:BALL] = bset = Set{Int64}()
for (elid, eltype) in mesh.element_types
    if eltype == :Tet4
        push!(bset, elid)
    end
end

for (elset_name, element_ids) in mesh.element_sets
    nel = length(element_ids)
    println("Element set $elset_name contains $nel elements.")
end

for (nset_name, node_ids) in mesh.node_sets
    nno = length(node_ids)
    println("Node set $nset_name contains $nno nodes.")
end

nnodes = length(mesh.nodes)
println("Total number of nodes in mesh: $nnodes")
nelements = length(mesh.elements)
println("Total number of elements in mesh: $nelements")

ball_elements = create_elements(mesh, "BALL")
nel_ball_elements = length(ball_elements)
println("ball contains $nel_ball_elements elements")
update!(ball_elements, "youngs modulus", 288.0)
update!(ball_elements, "poissons ratio", 1/3)
update!(ball_elements, "density", 36.0e-3)
#update!(ball_elements, "displacement load 1", 0.20)
update!(ball_elements, "displacement load 2", 0.20)

xmax = -Inf
nid = 0
for (j, coord) in mesh.nodes
    x, y, z = coord
    if x > xmax
        xmax = x
        nid = j
    end
end
println("max_x node = $nid, xval = $xmax")

load_element = Element(Poi1, [nid])
update!(load_element, "geometry", mesh.nodes)

V_ball = 0.0
m_ball = 0.0
time = 0.0
for element in ball_elements
    for ip in get_integration_points(element)
        detJ = element(ip, time, Val{:detJ})
        rho = element("density", ip, time)
        V_ball += ip.weight * detJ
        m_ball += ip.weight * rho * detJ
    end
end
println("Ball volume: $V_ball")
println("Ball mass: $m_ball")

ball = Problem(Elasticity, "ball", 3)
add_elements!(ball, ball_elements)
update!(load_element, "displacement traction force z", 10e4)

# Model construction end.

analysis = Analysis(LinearImplicit)
analysis.properties.tspan = (0.0, 10.0)
add_problems!(analysis, [ball])
xdmf = Xdmf(joinpath(datadir, "ball_results"); overwrite=true)
add_results_writer!(analysis, xdmf)

# Start analysis.

run!(analysis)

# Analysis ready

# sol = analysis.properties.sol
#
# dim = 3
# nnodes = length(mesh.nodes)
# ndofs = dim * nnodes
#
# u = hcat([x[1:ndofs] for x in sol.u]...)
# v = hcat([x[ndofs+1:end] for x in sol.u]...)
#
# tspan = analysis.properties.tspan
# t0, t1 = tspan
# nsteps = length(sol.u)
# println("Number of time steps = $nsteps")
# t = linspace(t0, t1, nsteps)
#
# for (time, x) in zip(t, sol.u)
#     u = x[1:ndofs]
#     v = x[ndofs+1:end]
#     ur = reshape(u, dim, nnodes)
#     vr = reshape(v, dim, nnodes)
#     ud = Dict(j => ur[:,j] for j=1:nnodes)
#     vd = Dict(j => vr[:,j] for j=1:nnodes)
#     update!(ball_elements, "displacement", time => ud)
#     update!(ball_elements, "velocity", time => vd)
# end

# nid2 = find_nearest_node(mesh, [0.0, 0.0, 0.0])

# ball("displacement", 0.0)
# trajectory = [ball("displacement", time)[nid2] for time in t]
# traj_x = [t[1] for t in trajectory]
# traj_y = [t[2] for t in trajectory]
# traj_z = [t[3] for t in trajectory]
# xmin, xmax = minimum(traj_x), maximum(traj_x)
# ymin, ymax = minimum(traj_y), maximum(traj_y)
# zmin, zmax = minimum(traj_z), maximum(traj_z)
# println("bounding box x = ($xmin, $xmax)")
# println("bounding box y = ($ymin, $ymax)")
# println("bounding box z = ($zmin, $zmax)")
#
# for time in linspace(0.0, 10.0, 600)
#     JuliaFEM.write_results!(analysis, time)
# end
close(xdmf.hdf) # src

#


text:
module Darknet
using Libdl, CEnum

using Darknet_jll

include(joinpath(@__DIR__, "..", "gen", "ctypes.jl"))
export Ctm, Ctime_t, Cclock_t
include(joinpath(@__DIR__, "..", "lib", "LibDarknet.jl"))

include("testing.jl")

include("utils.jl")
struct Image{N}
    w::Cint
    h::Cint
    c::Cint
    data::Array{Cfloat, N}
    Image(w,h,c,data::Array{Cfloat, N}) where N = new{N}(w,h,c,data)
end
Base.cconvert(::Type{image}, im::Image) = im
Base.unsafe_convert(::Type{image}, im::Image) = image(im.w, im.h, im.c, pointer(im.data))

end # module Darknet

ext: jl


text:
module TestElixirs

using LinearAlgebra
using Test
using Trixi

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
const EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples")

@testset "Special elixirs" begin
  @testset "Convergence test" begin
    mean_values = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "2d", "elixir_advection_extended.jl"), 3)
    @test isapprox(mean_values[:l2], [4.0], rtol=0.01)
  end

  @testset "Test linear structure (2D)" begin
    trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "2d", "elixir_advection_extended.jl"),
                  tspan=(0.0, 0.0), initial_refinement_level=2)
    A, b = linear_structure(semi)
    λ = eigvals(Matrix(A))
    @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))

    trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "2d", "elixir_hypdiff_lax_friedrichs.jl"),
                  tspan=(0.0, 0.0), initial_refinement_level=2)
    A, b = linear_structure(semi)
    λ = eigvals(Matrix(A))
    @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))
  end

  @testset "Test Jacobian of DG (2D)" begin
    trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "2d", "elixir_advection_extended.jl"),
                  tspan=(0.0, 0.0), initial_refinement_level=2)
    A, _ = linear_structure(semi)
    J = jacobian_fd(semi)
    @test Matrix(A) ≈ J
    λ = eigvals(J)
    @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))

    trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "2d", "elixir_euler_density_wave.jl"),
                  tspan=(0.0, 0.0), initial_refinement_level=2)
    J = jacobian_fd(semi)
    λ = eigvals(J)
    @test maximum(real, λ) < 0.007
  end

  @testset "Test linear structure (3D)" begin
    trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "3d", "elixir_advection_extended.jl"),
                  tspan=(0.0, 0.0), initial_refinement_level=1)
    A, b = linear_structure(semi)
    λ = eigvals(Matrix(A))
    @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))
  end

  @testset "Test Jacobian of DG (3D)" begin
    trixi_include(@__MODULE__, joinpath(EXAMPLES_DIR, "3d", "elixir_advection_extended.jl"),
                  tspan=(0.0, 0.0), initial_refinement_level=1)
    A, _ = linear_structure(semi)
    J = jacobian_fd(semi)
    @test Matrix(A) ≈ J
    λ = eigvals(J)
    @test maximum(real, λ) < 10 * sqrt(eps(real(semi)))
  end
end

# Clean up afterwards: delete Trixi output directory
@test_nowarn rm(outdir, recursive=true)

end #module


text:
"""
Normal Mixture of `k` components via an online EM algorithm.  `start` is a keyword
argument specifying the initial parameters.

```julia
o = NormalMix(2, LearningRate(); start = MixtureModel(Normal, [(0, 1), (3, 1)]))
mean(o)
var(o)
std(o)
```
"""
type NormalMix{W<:Weight} <: DistributionStat{ScalarInput}
    value::Ds.MixtureModel{Ds.Univariate, Ds.Continuous, Ds.Normal{Float64}}
    s1::VecF
    s2::VecF
    s3::VecF
    w::VecF
    μ::VecF
    σ2::VecF
    weight::W
end
function NormalMix(k::Integer, wgt::Weight = LearningRate();
        start = Ds.MixtureModel(Ds.Normal{Float64}[Ds.Normal(j-1, 10) for j in 1:k])
    )
    NormalMix(
        start,
        ones(k) / k, zeros(k), zeros(k),  # s1, s2, s3
        zeros(k), collect(1.:k), fill(10.0, k),  # w, μ, σ2
        wgt
    )
end
function NormalMix(y::AVec, k::Integer, wgt::Weight = LearningRate(); kw...)
    o = NormalMix(k, wgt; kw...)
    fit!(o, y)
    o
end
NormalMix(k::Integer, y::AVec, wgt::Weight = LearningRate(); kw...) = NormalMix(y, k, wgt; kw...)
Ds.componentwise_pdf(o::NormalMix, y) = Ds.componentwise_pdf(value(o), y)
Ds.ncomponents(o::NormalMix) = Ds.ncomponents(value(o))
Ds.component(o::NormalMix, j) = Ds.component(value(o), j)
Ds.probs(o::NormalMix) = Ds.probs(value(o))
Ds.pdf(o::NormalMix, y) = Ds.pdf(value(o), y)
Ds.cdf(o::NormalMix, y) = Ds.cdf(value(o), y)
Base.mean(o::NormalMix) = mean(value(o))
Base.var(o::NormalMix) = var(value(o))
Base.std(o::NormalMix) = std(value(o))
function Base.show(io::IO, o::NormalMix)
    printheader(io, "NormalMix (k = $(Ds.ncomponents(o)))")
    print_value_and_nobs(io, o)
end
function value(o::NormalMix)
    # try
        vec = map((u,v) -> Ds.Normal{Float64}(u, sqrt(v)), o.μ, o.σ2)
        o.value = Ds.MixtureModel(vec, o.s1)
    # catch
    #     println(sqrt(o.σ2))
    #     println(o.μ)
    #     println(o.s1)
    #     error("Algorithm possibly diverging, nobs = $(nobs(o))")
    # end
end
function _fit!(o::NormalMix, y::Real, γ::Float64)
    k = length(o.μ)
    for j in 1:k
        σinv = 1.0 / sqrt(o.σ2[j])
        o.w[j] = o.s1[j] * σinv * exp(-.5 * σinv * σinv * (y - o.μ[j]) ^ 2)
    end
    sum1 = sum(o.w)
    for j in 1:k
        o.w[j] /= sum1
        o.s1[j] = smooth(o.s1[j], o.w[j], γ)
        o.s2[j] = smooth(o.s2[j], o.w[j] * y, γ)
        o.s3[j] = smooth(o.s3[j], o.w[j] * y * y, γ)
    end
    sum2 = sum(o.s1)
    for j in 1:k
        o.μ[j] = o.s2[j] / o.s1[j]
        o.σ2[j] = (o.s3[j] - o.s2[j] ^ 2 / o.s1[j]) / o.s1[j]
        o.s1[j] /= sum2
        if o.σ2[j] <= _ϵ
            o.σ2 = ones(k)
        end
    end
    o
end
function _fitbatch!{T<:Real}(o::NormalMix, y::AVec{T}, γ::Float64)
    n = length(y)
    k = length(o.μ)
    s1 = copy(o.s1)
    s2 = copy(o.s2)
    s3 = copy(o.s3)

    for yi in y
        for j in 1:k
            σinv = 1.0 / sqrt(o.σ2[j])
            o.w[j] = s1[j] * σinv * exp(-.5 * σinv * σinv * (yi - o.μ[j]) ^ 2)
        end
        sum1 = sum(o.w)
        for j in 1:k
            o.w[j] /= sum1
            if yi == y[1]
                o.s1[j] = smooth(o.s1[j], o.w[j] / n, γ)
                o.s2[j] = smooth(o.s2[j], o.w[j] * yi / n, γ)
                o.s3[j] = smooth(o.s3[j], o.w[j] * yi * yi / n, γ)
            else
                o.s1[j] += γ * o.w[j] / n
                o.s2[j] += γ * o.w[j] * yi / n
                o.s3[j] += γ * o.w[j] * yi * yi / n
            end
        end
    end
    sum2 = sum(o.s1)
    for j in 1:k
        o.μ[j] = o.s2[j] / o.s1[j]
        o.σ2[j] = (o.s3[j] - o.s2[j] ^ 2 / o.s1[j]) / o.s1[j]
        o.s1[j] /= sum2
        if o.σ2[j] <= _ϵ
            o.σ2 = ones(k)
        end
    end
    o
end

# Quantiles.  Starting values based on Normal distribution.
function Base.quantile(o::NormalMix, τ::Real;
        start = quantile(Ds.Normal(mean(o), std(o)), τ), maxit = 20, tol = .001
    )
    @assert 0 < τ < 1
    θ = start
    for i in 1:maxit
        θ += (τ - Ds.cdf(o, θ)) / Ds.pdf(o, θ)
        abs(Ds.cdf(o, θ) - τ) < tol && break
    end
    return θ
end
function Base.quantile{T<:Real}(o::NormalMix, τ::Vector{T};
        start = start = quantile(Ds.Normal(mean(o), std(o)), τ), kw...
    )
    Float64[quantile(o, τ[j]; start = start[j], kw...) for j in 1:length(τ)]
end

# #-----------------------------------------------------------------# Multivariate
# type MvNormalMix{W<:Weight} <: OnlineStat
#     value::Ds.MixtureModel{Ds.Multivariate, Ds.Continuous, Ds.MvNormal}
#     s1::VecF
#     s2::MatF
#     s3::Array{Float64, 3}
#     w::VecF
#     μ::MatF
#     V::Array{Float64, 3}
#     weight::W
#     n::Int
#     nups::Int
# end
# function MvNormalMix(p::Integer, k::Integer, wgt::Weight = LearningRate();
#         start = Ds.MixtureModel(Ds.MvNormal[Ds.MvNormal(ones(p) + j, 10 * eye(p)) for j in 1:k])
#     )
#     V = zeros(p, p, k)
#     # for j in 1:k
#     #     V[:, :, j] = cov(Ds.component(start, j))
#     # end
#
#     MvNormalMix(start,
#         ones(k) / k, zeros(p, k), zeros(p, p, k),  # s1, s2, s3
#         zeros(k), zeros(p, k), V,            # w, μ, V
#         wgt, 0, 0
#     )
# end
# function MvNormalMix{T<:Real}(y::AMat{T}, k::Integer, wgt::Weight = LearningRate(); kw...)
#     o = MvNormalMix(size(y, 2), k, wgt; kw...)
#     fit!(o, y)
#     o
# end
#
# function fit!{T<:Real}(o::MvNormalMix, y::AVec{T})
#     k = length(o.w)
#     γ = weight!(o, 1)
#     for j in 1:k
#         o.w[j] = o.s1[j] * Ds.pdf(Ds.component(o.value, j), y)
#     end
#     sum1 = sum(o.w)
#     for j in 1:k
#         o.w[j] /= sum1
#         o.s1[j] = smooth(o.s1[j], o.w[j], γ)
#         smooth!(o.s2[:, j], o.w[j] * y, γ)
#         smooth!(o.s3[:, :, j], o.w[j] * y * y', γ)
#
#         o.μ[:, j] = o.s2[j] / o.s1[j]
#         o.V[:, :, j] = (o.s3[:, :, j] - o.s2[j] * o.s2[j]' / o.s1[j]) / o.s1[j]
#     end
#
#     sum2 = sum(o.s1)
#     for j in 1:k
#         o.s1[j] /= sum2
#     end
#
#     o.value = Ds.MixtureModel(
#         Ds.MvNormal[
#             Ds.MvNormal(o.μ[:,j], o.V[:,:,j] + .001 * eye(length(o.μ[:,j]))) for j in 1:k
#         ],
#         o.s1 / sum(o.s1)
#     )
# end
#


text:
@testset "Peak width" begin
    fs = 100
    T = 1/fs
    t = 0:T:6pi+T
    sint = sin.(t)

    @testset "minima/maxima" begin
        sinpks = argmaxima(sint)
        _, widths, _, _ = peakwidths(sinpks, sint, sint[sinpks]; strict=false, relheight=1)
        @test widths ≈ fill(pi*100, length(sinpks)) atol=.01

        sinpks = argminima(sint)
        _, widths, _, _ = peakwidths(sinpks, sint, abs.(sint[sinpks]); strict=false, relheight=1)
        @test widths ≈ fill(pi*100, length(sinpks)) atol=.01
    end

    _, widths, _, _ = peakwidths([2], [0.,1.,0.], [1.])
    @test widths == [1.]
    _, widths, _, _ = peakwidths([2], [0.,1.,0.], [NaN])
    @test widths[1] === NaN
    _, widths, _, _ = peakwidths([2], [0.,1.,0.], [missing])
    @test widths[1] === missing
    _, widths, _, _ = peakwidths([2], [0.,1.,NaN], [1.]; strict=true)
    @test widths[1] === NaN
    _, widths, _, _ = peakwidths([2], [0.,1.,0.,-1.], [1.]; strict=false)
    _, widthsnan, _, _ = peakwidths([2], [0.,1.,NaN,-1.], [1.]; strict=false)
    @test widths == widthsnan

    @testset "strict" begin
        sinpks = argmaxima(sint)
        _, widths, _, _ = peakwidths(sinpks, sint, ones(length(sinpks)); strict=true, relheight=1)
        @test first(widths) === NaN
        _, widths, _, _ = peakwidths(sinpks, sint, ones(length(sinpks)); strict=false, relheight=1)
        @test first(widths) !== NaN
    end

    @testset "Min/max width" begin
        sinpks = argmaxima(sint)
        _, proms = peakproms(sinpks, sint)

        @test length(first(peakwidths(sinpks, sint, proms; minwidth=pi*75))) == 2
        @test length(first(peakwidths(sinpks, sint, proms; maxwidth=pi*75))) == 1

        @test_throws ArgumentError peakwidths(1:3, ones(3), ones(3); maxwidth=0.1, minwidth=1)
    end

end

max_stars_repo_name: halleysfifthinc/Peaks

max_issues_repo_name: halleysfifthinc/Peaks

max_forks_repo_name: halleysfifthinc/Peaks


text:
using .NESInterface, Cairo

include("actions.jl")
include("decode_target.jl")

mutable struct SMBEnv

    target_world
    target_area
    target_stage

    time_last
    x_position_last
    nesenv::NESInterface.NESEnv

    action_map::Dict{Int, UInt8}
    action_meanings::Dict{Int, String}
end

const SMB_reward_range = (-15, 15)

function SMBEnv(rom_path::AbstractString, action_type::Symbol, lost_levels=false, target=nothing)

    target_world, target_stage, target_area = decode_target(target, lost_levels)

    actions = action_types_dict[action_type]
    time_last = 0
    x_position_last = 0

    nesenv = NESInterface.NESEnv(rom_path)

    action_map = Dict{Int, UInt8}()
    action_meanings = Dict{Int, String}()

    for (action, button_list) ∈ enumerate(actions)
        byte_action = 0
        for  button ∈ button_list
            byte_action |= button_map[button]
        end
        action_map[action] = byte_action
        action_meanings[action] = join(button_list, ' ')
    end

    SMBEnv(target_world, target_area, target_stage, time_last, x_position_last,
            nesenv, action_map, action_meanings)
end

"""Read a range of bytes where each byte is a 10's place figure.
Note:
    this method is specific to Mario where three GUI values are stored
    in independent memory slots to save processing time
    - score has 6 10's places
    - coins has 2 10's places
    - time has 3 10's places
"""
function read_mem_range(smb::SMBEnv, address, length)
    accum = 0
    for i=address:address+length-1
        accum *= 10
        accum += smb.nesenv.ram[i]
    end
    return accum
end

function Base.getproperty(smb::SMBEnv, sym::Symbol)
    if sym == :is_single_stage_env
        return !isnothing(smb.target_world) && !isnothing(smb.target_area)

    elseif sym == :level
        return smb.nesenv.ram[0x7560] * 4 + smb.nesenv.ram[0x075d]

    elseif sym == :world
        return smb.nesenv.ram[0x0760] + 1

    elseif sym == :stage
        return smb.nesenv.ram[0x075d] + 1

    elseif sym == :area
        return smb.nesenv.ram[0x0761] + 1

    elseif sym == :score
        return read_mem_range(smb, 0x07dd, 6)

    elseif sym == :time
        return read_mem_range(smb, 0x07f9, 3)

    elseif sym == :coins
        return read_mem_range(smb, 0x07ee, 2)

    elseif sym == :life
        return smb.nesenv.ram[0x075b]

    elseif sym == :x_position
        return smb.nesenv.ram[0x6e] * 0x100 + smb.nesenv.ram[0x87]

    elseif sym == :left_x_position
        return (smb.nesenv.ram[0x87] - smb.nesenv.ram[0x071d]) % 256

    elseif sym == :y_position
        return smb.nesenv.ram[0x03b9]

    elseif sym == :y_viewport
        return smb.nesenv.ram[0x00b6]

    elseif sym == :player_status
        status = smb.nesenv.ram[0x0757]

        if status == 0
            return :small
        elseif status == 1
            return :tall
        elseif status == 2
            return :fireball
        end

    elseif sym == :player_state
        return smb.nesenv.ram[0x000f]

    elseif sym == :is_dying
        return smb.player_state == 0x0b || smb.y_viewport > 1

    elseif sym == :is_dead
        return smb.player_state == 0x06

    elseif sym == :is_game_over
        return smb.life == 0xff

    elseif sym == :is_busy
        return smb.player_state ∈ (0x00, 0x01, 0x02, 0x03, 0x04, 0x05, 0x07)

    elseif sym == :is_world_over
        return smb.nesenv.ram[0x0771] == 2

    elseif sym == :is_stage_over
        for address ∈ (0x0017, 0x0018, 0x0019, 0x001a, 0x001b)
            if smb.nesenv.ram[address] ∈ (0x2d, 0x31)
                return smb.nesenv.ram[0x001e] == 3
            end
        end
        return false

    elseif sym == :flag_get
        return smb.is_world_over || smb.is_stage_over

    elseif sym == :death_penalty
        return smb.is_dying || smb.is_dead ? -25 : 0

    elseif sym == :x_reward
        _reward = smb.x_position - smb.x_position_last
        smb.x_position_last = smb.x_position

        return -5 ≤ _reward ≤ 5 ? _reward : 0

    elseif sym == :time_penalty
        _reward = smb.time - smb.time_last
        smb.time_last = smb.time
        # time can only decrease, a positive reward results from a reset an
        # should default to 0 reward
        return min(0, _reward)

    else
        return Base.getfield(smb, sym)
    end
end

"""Force the prelevel timer to 0 to skip frams during a death"""
runout_prelevel_timer!(smb::SMBEnv) = (smb.nesenv.ram[0x07a1] = 0)

# ------------------------ Reset level and its utilites ------------------------
"""Write the stage data to RAM to overwrite loading the next stage."""
function write_stage!(smb::SMBEnv)
    smb.nesenv.ram[0x0760] = smb.target_world - 1
    smb.nesenv.ram[0x075d] = smb.target_stage - 1
    smb.nesenv.ram[0x0761] = smb.target_area - 1
end

"""Press and release start to skip the start screen"""
function skip_start_screen!(smb::SMBEnv)
    # Press and release the start button
    NESInterface.frame_advance!(smb.nesenv, 8)
    NESInterface.frame_advance!(smb.nesenv, 0)

    # Press start until the game starts
    while smb.time == 0
        # Press and release the start button
        NESInterface.frame_advance!(smb.nesenv, 8)
        # if we're in the single stage environment, write the stage data
        if smb.is_single_stage_env
            write_stage!(smb)
        end
        NESInterface.frame_advance!(smb.nesenv, 0)
        # run-out the prelevel times to skip the animation
        runout_prelevel_timer!(smb)
    end

    # set the last time to now
    smb.time_last = smb.time
    # after the start screen idle to skip some extra frames
    while smb.time >= smb.time_last
        smb.time_last = smb.time
        NESInterface.frame_advance!(smb.nesenv, 8)
        NESInterface.frame_advance!(smb.nesenv, 0)
    end
end

function reset!(smb::SMBEnv)
    smb.time_last = 0
    smb.x_position_last = 0
    NESInterface.reset!(smb.nesenv)
    skip_start_screen!(smb)
    smb.time_last = smb.time
    smb.x_position_last = smb.x_position
    smb.nesenv.screen
end

# --------------------- Step Function and its utils ----------------------------

get_reward(smb::SMBEnv) = smb.x_reward + smb.time_penalty + smb.death_penalty
get_done(smb::SMBEnv) = smb.is_single_stage_env ? smb.is_dying || smb.is_dead || smb.flag_get :
                            smb.is_game_over
function get_info(smb::SMBEnv)
    return Dict(
        :coins    => smb.coins,
        :flag_get => smb.flag_get,
        :life     => Int(smb.life),
        :score    => smb.score,
        :stage    => smb.stage,
        :status   => smb.player_status,
        :time     => smb.time,
        :world    => smb.world,
        :x_pos    => Int(smb.x_position)
    )
end

"""If Mario dies, then skip the death animation..."""
function kill_mario!(smb::SMBEnv)
    smb.nesenv.ram[0x00f] = 0x06
    NESInterface.frame_advance!(smb.nesenv, 0)
end

"""Skip the cutscene that plays at the end of a world"""
function skip_end_of_world!(smb::SMBEnv)
    if smb.is_world_over
        # Get the current game time to reference.
        time = smb.time
        # Loop until the time is different
        while smb.time == time
            NESInterface.frame_advance!(smb.nesenv, 0)
        end
    end
end

"""Skip change area animations by running down timers"""
function skip_area_change!(smb::SMBEnv)
    change_area_timer = smb.nesenv.ram[0x006df]
    1 < change_area_timer < 255 && (smb.nesenv.ram[0x06df] = 1)
end

"""Skip occupied states by running out a timer and skipping frames."""
function skip_occupied_states!(smb::SMBEnv)
    while smb.is_busy || smb.is_world_over
        runout_prelevel_timer!(smb)
        NESInterface.frame_advance!(smb.nesenv, 0)
    end
end

function step!(smb::SMBEnv, action)
    smb.nesenv.done && error("Cannot step in a completed environment! Call `reset!`")
    if action ∉ keys(smb.action_map)
        error("Invalid action entered.")
    end

    # Perform step
    NESInterface.frame_advance!(smb.nesenv, smb.action_map[action])

    # Calculate reward, done and step info
    reward = get_reward(smb)
    smb.nesenv.done = get_done(smb)
    info = get_info(smb)

    # Post step hacking
    if !smb.nesenv.done
        smb.is_dying && kill_mario!(smb)
        smb.is_single_stage_env || skip_end_of_world!(smb)

        # skip area change (i.e. enter pipe, flag get, etc.)
        skip_area_change!(smb)
        # skip occupied states like the black screen between lives that shows how
        # many lives the player has left
        skip_occupied_states!(smb)
    end

    reward = clamp(reward, SMB_reward_range[1], SMB_reward_range[2])
    return smb.nesenv.screen, reward, smb.nesenv.done, info
end

function render(smb::SMBEnv)
    screen = smb.nesenv.screen
    uint32_arr = screen[:, :, 1] .+ screen[:, :, 2] * 0x00000100 .+ screen[:, :, 3] * 0x00010000
    uint32_arr = uint32_arr |> transpose |> Array
    CairoRGBSurface(uint32_arr)
end


text:
# Automatically-generated headers for nvtx
#
# DO NOT EDIT THIS FILE DIRECTLY.

# Julia wrapper for header: nvToolsExt.h
# Automatically generated using Clang.jl

function nvtxInitialize(reserved)
    initialize_api()
    @runtime_ccall((:nvtxInitialize, libnvtx()), Cvoid,
                   (Ptr{Cvoid},),
                   reserved)
end

function nvtxDomainMarkEx(domain, eventAttrib)
    initialize_api()
    @runtime_ccall((:nvtxDomainMarkEx, libnvtx()), Cvoid,
                   (nvtxDomainHandle_t, Ptr{nvtxEventAttributes_t}),
                   domain, eventAttrib)
end

function nvtxMarkEx(eventAttrib)
    initialize_api()
    @runtime_ccall((:nvtxMarkEx, libnvtx()), Cvoid,
                   (Ptr{nvtxEventAttributes_t},),
                   eventAttrib)
end

function nvtxMarkA(message)
    initialize_api()
    @runtime_ccall((:nvtxMarkA, libnvtx()), Cvoid,
                   (Cstring,),
                   message)
end

function nvtxMarkW(message)
    initialize_api()
    @runtime_ccall((:nvtxMarkW, libnvtx()), Cvoid,
                   (Ptr{Cwchar_t},),
                   message)
end

function nvtxDomainRangeStartEx(domain, eventAttrib)
    initialize_api()
    @runtime_ccall((:nvtxDomainRangeStartEx, libnvtx()), nvtxRangeId_t,
                   (nvtxDomainHandle_t, Ptr{nvtxEventAttributes_t}),
                   domain, eventAttrib)
end

function nvtxRangeStartEx(eventAttrib)
    initialize_api()
    @runtime_ccall((:nvtxRangeStartEx, libnvtx()), nvtxRangeId_t,
                   (Ptr{nvtxEventAttributes_t},),
                   eventAttrib)
end

function nvtxRangeStartA(message)
    initialize_api()
    @runtime_ccall((:nvtxRangeStartA, libnvtx()), nvtxRangeId_t,
                   (Cstring,),
                   message)
end

function nvtxRangeStartW(message)
    initialize_api()
    @runtime_ccall((:nvtxRangeStartW, libnvtx()), nvtxRangeId_t,
                   (Ptr{Cwchar_t},),
                   message)
end

function nvtxDomainRangeEnd(domain, id)
    initialize_api()
    @runtime_ccall((:nvtxDomainRangeEnd, libnvtx()), Cvoid,
                   (nvtxDomainHandle_t, nvtxRangeId_t),
                   domain, id)
end

function nvtxRangeEnd(id)
    initialize_api()
    @runtime_ccall((:nvtxRangeEnd, libnvtx()), Cvoid,
                   (nvtxRangeId_t,),
                   id)
end

function nvtxDomainRangePushEx(domain, eventAttrib)
    initialize_api()
    @runtime_ccall((:nvtxDomainRangePushEx, libnvtx()), Cint,
                   (nvtxDomainHandle_t, Ptr{nvtxEventAttributes_t}),
                   domain, eventAttrib)
end

function nvtxRangePushEx(eventAttrib)
    initialize_api()
    @runtime_ccall((:nvtxRangePushEx, libnvtx()), Cint,
                   (Ptr{nvtxEventAttributes_t},),
                   eventAttrib)
end

function nvtxRangePushA(message)
    initialize_api()
    @runtime_ccall((:nvtxRangePushA, libnvtx()), Cint,
                   (Cstring,),
                   message)
end

function nvtxRangePushW(message)
    initialize_api()
    @runtime_ccall((:nvtxRangePushW, libnvtx()), Cint,
                   (Ptr{Cwchar_t},),
                   message)
end

function nvtxDomainRangePop(domain)
    initialize_api()
    @runtime_ccall((:nvtxDomainRangePop, libnvtx()), Cint,
                   (nvtxDomainHandle_t,),
                   domain)
end

function nvtxRangePop()
    initialize_api()
    @runtime_ccall((:nvtxRangePop, libnvtx()), Cint, ())
end

function nvtxDomainResourceCreate(domain, attribs)
    initialize_api()
    @runtime_ccall((:nvtxDomainResourceCreate, libnvtx()), nvtxResourceHandle_t,
                   (nvtxDomainHandle_t, Ptr{nvtxResourceAttributes_t}),
                   domain, attribs)
end

function nvtxDomainResourceDestroy(resource)
    initialize_api()
    @runtime_ccall((:nvtxDomainResourceDestroy, libnvtx()), Cvoid,
                   (nvtxResourceHandle_t,),
                   resource)
end

function nvtxDomainNameCategoryA(domain, category, name)
    initialize_api()
    @runtime_ccall((:nvtxDomainNameCategoryA, libnvtx()), Cvoid,
                   (nvtxDomainHandle_t, UInt32, Cstring),
                   domain, category, name)
end

function nvtxDomainNameCategoryW(domain, category, name)
    initialize_api()
    @runtime_ccall((:nvtxDomainNameCategoryW, libnvtx()), Cvoid,
                   (nvtxDomainHandle_t, UInt32, Ptr{Cwchar_t}),
                   domain, category, name)
end

function nvtxNameCategoryA(category, name)
    initialize_api()
    @runtime_ccall((:nvtxNameCategoryA, libnvtx()), Cvoid,
                   (UInt32, Cstring),
                   category, name)
end

function nvtxNameCategoryW(category, name)
    initialize_api()
    @runtime_ccall((:nvtxNameCategoryW, libnvtx()), Cvoid,
                   (UInt32, Ptr{Cwchar_t}),
                   category, name)
end

function nvtxNameOsThreadA(threadId, name)
    initialize_api()
    @runtime_ccall((:nvtxNameOsThreadA, libnvtx()), Cvoid,
                   (UInt32, Cstring),
                   threadId, name)
end

function nvtxNameOsThreadW(threadId, name)
    initialize_api()
    @runtime_ccall((:nvtxNameOsThreadW, libnvtx()), Cvoid,
                   (UInt32, Ptr{Cwchar_t}),
                   threadId, name)
end

function nvtxDomainRegisterStringA(domain, string)
    initialize_api()
    @runtime_ccall((:nvtxDomainRegisterStringA, libnvtx()), nvtxStringHandle_t,
                   (nvtxDomainHandle_t, Cstring),
                   domain, string)
end

function nvtxDomainRegisterStringW(domain, string)
    initialize_api()
    @runtime_ccall((:nvtxDomainRegisterStringW, libnvtx()), nvtxStringHandle_t,
                   (nvtxDomainHandle_t, Ptr{Cwchar_t}),
                   domain, string)
end

function nvtxDomainCreateA(name)
    initialize_api()
    @runtime_ccall((:nvtxDomainCreateA, libnvtx()), nvtxDomainHandle_t,
                   (Cstring,),
                   name)
end

function nvtxDomainCreateW(name)
    initialize_api()
    @runtime_ccall((:nvtxDomainCreateW, libnvtx()), nvtxDomainHandle_t,
                   (Ptr{Cwchar_t},),
                   name)
end

function nvtxDomainDestroy(domain)
    initialize_api()
    @runtime_ccall((:nvtxDomainDestroy, libnvtx()), Cvoid,
                   (nvtxDomainHandle_t,),
                   domain)
end
# Julia wrapper for header: nvToolsExtCuda.h
# Automatically generated using Clang.jl

function nvtxNameCuDeviceA(device, name)
    initialize_api()
    @runtime_ccall((:nvtxNameCuDeviceA, libnvtx()), Cvoid,
                   (CUdevice, Cstring),
                   device, name)
end

function nvtxNameCuDeviceW(device, name)
    initialize_api()
    @runtime_ccall((:nvtxNameCuDeviceW, libnvtx()), Cvoid,
                   (CUdevice, Ptr{Cwchar_t}),
                   device, name)
end

function nvtxNameCuContextA(context, name)
    initialize_api()
    @runtime_ccall((:nvtxNameCuContextA, libnvtx()), Cvoid,
                   (CUcontext, Cstring),
                   context, name)
end

function nvtxNameCuContextW(context, name)
    initialize_api()
    @runtime_ccall((:nvtxNameCuContextW, libnvtx()), Cvoid,
                   (CUcontext, Ptr{Cwchar_t}),
                   context, name)
end

function nvtxNameCuStreamA(stream, name)
    initialize_api()
    @runtime_ccall((:nvtxNameCuStreamA, libnvtx()), Cvoid,
                   (CUstream, Cstring),
                   stream, name)
end

function nvtxNameCuStreamW(stream, name)
    initialize_api()
    @runtime_ccall((:nvtxNameCuStreamW, libnvtx()), Cvoid,
                   (CUstream, Ptr{Cwchar_t}),
                   stream, name)
end

function nvtxNameCuEventA(event, name)
    initialize_api()
    @runtime_ccall((:nvtxNameCuEventA, libnvtx()), Cvoid,
                   (CUevent, Cstring),
                   event, name)
end

function nvtxNameCuEventW(event, name)
    initialize_api()
    @runtime_ccall((:nvtxNameCuEventW, libnvtx()), Cvoid,
                   (CUevent, Ptr{Cwchar_t}),
                   event, name)
end


text:
using ITensorGaussianMPS
using ITensors
using LinearAlgebra
using Test

@testset "Basic" begin
  # Test Givens rotations
  v = randn(6)
  g, r = ITensorGaussianMPS.givens_rotations(v)
  @test g * v ≈ r * [n == 1 ? 1 : 0 for n in 1:length(v)]
end

@testset "Fermion" begin
  N = 10
  Nf = N ÷ 2

  # Hopping
  t = 1.0

  # Hopping Hamiltonian
  h = Hermitian(diagm(1 => fill(-t, N - 1), -1 => fill(-t, N - 1)))
  e, u = eigen(h)

  @test h * u ≈ u * Diagonal(e)

  E = sum(e[1:Nf])

  # Get the Slater determinant
  Φ = u[:, 1:Nf]
  @test h * Φ ≈ Φ * Diagonal(e[1:Nf])

  # Diagonalize the correlation matrix as a
  # Gaussian MPS (GMPS) gates
  n, gmps = slater_determinant_to_gmera(Φ; maxblocksize=10)

  ns = round.(Int, n)
  @test sum(ns) == Nf

  Λ = conj(Φ) * transpose(Φ)
  @test gmps * Λ * gmps' ≈ Diagonal(ns) rtol = 1e-2
  @test gmps' * Diagonal(ns) * gmps ≈ Λ rtol = 1e-2

  # Form the MPS
  s = siteinds("Fermion", N; conserve_qns=true)
  ψ = ITensorGaussianMPS.slater_determinant_to_mera(s, Φ; blocksize=4)

  os = OpSum()
  for i in 1:N, j in 1:N
    if h[i, j] ≠ 0
      os .+= h[i, j], "Cdag", i, "C", j
    end
  end
  H = MPO(os, s)

  @test inner(ψ, H, ψ) ≈ E rtol = 1e-5

  # Compare to DMRG
  sweeps = Sweeps(10)
  setmaxdim!(sweeps, 10, 20, 40, 60)
  setcutoff!(sweeps, 1E-12)
  energy, ψ̃ = dmrg(H, productMPS(s, n -> n ≤ Nf ? "1" : "0"), sweeps; outputlevel=0)

  # Create an mps
  @test abs(inner(ψ, ψ̃)) ≈ 1 rtol = 1e-5
  @test inner(ψ̃, H, ψ̃) ≈ inner(ψ, H, ψ) rtol = 1e-5
  @test E ≈ energy
end

@testset "Fermion (complex)" begin
  N = 10
  Nf = N ÷ 2

  # Hopping
  θ = π / 8
  t = exp(im * θ)

  # Hopping Hamiltonian
  h = Hermitian(diagm(1 => fill(-t, N - 1), -1 => fill(-conj(t), N - 1)))
  e, u = eigen(h)

  @test h * u ≈ u * Diagonal(e)

  E = sum(e[1:Nf])

  # Get the Slater determinant
  Φ = u[:, 1:Nf]
  @test h * Φ ≈ Φ * Diagonal(e[1:Nf])

  # Diagonalize the correlation matrix as a
  # Gaussian MPS (GMPS)
  n, gmps = slater_determinant_to_gmera(Φ; maxblocksize=4)

  ns = round.(Int, n)
  @test sum(ns) == Nf

  Λ = conj(Φ) * transpose(Φ)
  @test gmps * Λ * gmps' ≈ Diagonal(ns) rtol = 1e-2
  @test gmps' * Diagonal(ns) * gmps ≈ Λ rtol = 1e-2

  # Form the MPS
  s = siteinds("Fermion", N; conserve_qns=true)
  ψ = ITensorGaussianMPS.slater_determinant_to_mera(s, Φ; blocksize=4)

  os = OpSum()
  for i in 1:N, j in 1:N
    if h[i, j] ≠ 0
      os .+= h[i, j], "Cdag", i, "C", j
    end
  end
  H = MPO(os, s)

  @test inner(ψ, H, ψ) ≈ E rtol = 1e-5
  @test inner(ψ, H, ψ) / norm(ψ) ≈ E rtol = 1e-5

  # Compare to DMRG
  sweeps = Sweeps(10)
  setmaxdim!(sweeps, 10, 20, 40, 60)
  setcutoff!(sweeps, 1E-12)
  energy, ψ̃ = dmrg(H, productMPS(s, n -> n ≤ Nf ? "1" : "0"), sweeps; outputlevel=0)

  # Create an mps
  @test abs(inner(ψ, ψ̃)) ≈ 1 rtol = 1e-5
  @test inner(ψ̃, H, ψ̃) ≈ inner(ψ, H, ψ) rtol = 1e-5
  @test E ≈ energy
end

# Build 1-d SSH model
function SSH1dModel(N::Int, t::Float64, vardelta::Float64)
  # N should be even
  s = siteinds("Fermion", N; conserve_qns=true)
  limit = div(N - 1, 2)
  t1 = -t * (1 + vardelta / 2)
  t2 = -t * (1 - vardelta / 2)
  os = OpSum()
  for n in 1:limit
    os .+= t1, "Cdag", 2 * n - 1, "C", 2 * n
    os .+= t1, "Cdag", 2 * n, "C", 2 * n - 1
    os .+= t2, "Cdag", 2 * n, "C", 2 * n + 1
    os .+= t2, "Cdag", 2 * n + 1, "C", 2 * n
  end
  if N % 2 == 0
    os .+= t1, "Cdag", N - 1, "C", N
    os .+= t1, "Cdag", N, "C", N - 1
  end
  h = hopping_hamiltonian(os)
  H = MPO(os, s)
  #display(t1)
  return (h, H, s)
end

@testset "Energy" begin
  N = 2^4
  Nf = div(N, 2)
  t = 1.0
  gapsize = 0
  vardelta = gapsize / 2
  h, H, s = SSH1dModel(N, t, vardelta)

  Φ = slater_determinant_matrix(h, Nf)
  E, V = eigen(h)
  sort(E)
  Eana = sum(E[1:Nf])

  Λ0 = Φ * Φ'
  @test Eana ≈ tr(h * Λ0) rtol = 1e-5
  # Diagonalize the correlation matrix as a
  # Gaussian MPS (GMPS) and GMERA
  ngmps, V1 = ITensorGaussianMPS.correlation_matrix_to_gmps(Λ0; eigval_cutoff=1e-8)
  nmera, V1 = ITensorGaussianMPS.correlation_matrix_to_gmera(Λ0; eigval_cutoff=1e-8)#,maxblocksize=6)
  @test sum(round.(Int, nmera)) == sum(round.(Int, ngmps))

  U = ITensorGaussianMPS.UmatFromGates(V1, N)
  Etest = ITensorGaussianMPS.EfromGates(h, U)

  @test Eana ≈ Etest rtol = 1e-5
end


text:
signp(x::Real) = x<0 ? -one(x) : one(x)

checkvec(n::Number) = n:n
checkvec(n::AbstractVector) = n

blasfloat(x::Int32) = Float32(x)
blasfloat(x::Int16) = Float32(x)
blasfloat(x::Number) = Float64(x)

matsize1(n,q) = 2 + maximum(n) + ceil(Int,sqrt(abs(q)))
matsize2(m,q) = 2 + maximum(m)÷2 + ceil(Int,sqrt(abs(q)))
matsize3(n,q) = 2 + maximum(n) + ceil(Int,2*sqrt(abs(q)))
matsize4(m,q) = 2 + maximum(m)÷2 + ceil(Int,2*sqrt(abs(q)))

Base.filter(::typeof(iseven), r::UnitRange{T}) where {T<:Integer} =
    StepRange{T,T}((first(r)+T(1)) & ~T(1), T(2), last(r) & ~T(1))
Base.filter(::typeof(isodd), r::UnitRange{T}) where {T<:Integer} =
    StepRange{T,T}(first(r) & ~T(1) + T(1), T(2), (last(r)+T(1)) & ~T(1) - T(1))


text:
import .Loggables

"""
    abstract type LoggerBackend

Backend for logging callbacks like.

To add support for logging [`Loggables.Loggable`](#) `L` to backend `B`, implement

[`log_to`](#)`(backend::B, loggable::L, names, i)`

See also [`LogMetrics`](#), [`LogHyperParams`](#), [`log_to`](#)
"""
abstract type LoggerBackend end

"""
    log_to(backend, loggable, group, i)
    log_to(backends, loggable, group, i)

Log `loggable` to `backend` with `group` to index `i`.

- `loggable` is any [`Loggables.Loggable`](#)
- `group` can be a `String` or a tuple of `String`s implying
  some grouping which can be used by a supporting backend.
- `i` is a step counter and unique for every group.
"""
function log_to(backends::Tuple, loggable, name, i; group = ())
    for backend in backends
        log_to(backend, loggable, name, i, group = group)
    end
end

"""
    LogMetrics(backends...) <: Callback

Callback that logs step and epoch metrics to one or more [`LoggerBackend`](#)s.

See also [`LoggerBackend`](#), [`Loggables.Loggable`](#), [`log_to`](#),
[`TensorBoardBackend`](#)

Example:

```julia
logcb = LogMetrics(TensorBoardBackend("tblogs"))
Learner(model, data, opt, lossfn, Metrics(accuracy), logcb)
```
"""
struct LogMetrics <: Callback
    backends::Tuple
    LogMetrics(backends...) = new(backends)
end

stateaccess(::LogMetrics) = (
    cbstate = (history = Read(), metricsstep = Read(), metricsepoch = Read()),)

function on(::StepEnd, phase, logger::LogMetrics, learner)
    history = learner.cbstate.history[phase]
    metricsstep = learner.cbstate.metricsstep[phase]
    for metric in keys(metricsstep)
        val = last(last(metricsstep, metric))
        log_to(
            logger.backends,
            Loggables.Value(val),
            string(metric),
            history.steps,
            group = ("Step", string(typeof(phase)), "Metrics"))
    end
end

function on(::EpochEnd, phase, logger::LogMetrics, learner)
    history = learner.cbstate.history[phase]
    metricsepoch = learner.cbstate.metricsepoch[phase]
    for metric in keys(metricsepoch)
        _, val = last(metricsepoch, metric)
        log_to(
            logger.backends,
            Loggables.Value(val),
            string(metric),
            history.epochs,
            group = ("Epoch", string(typeof(phase)), "Metrics"))
    end
end

"""
    LogHyperParams(backends...) <: Callback

Callback that logs hyperparameters to one or more [`LoggerBackend`](#)s.

See also [`LoggerBackend`](#), [`Loggables.Loggable`](#), [`log_to`](#),
[`TensorBoardBackend`](#)

## Example

```julia
logcb = LogHyperParams(TensorBoardBackend("tblogs"))
schedule = ...
Learner(model, data, opt, lossfn, Scheduler(LearningRate => schedule), logcb)
```
"""
struct LogHyperParams <: Callback
    backends::Tuple
    LogHyperParams(backends...) = new(backends)
end

stateaccess(::LogHyperParams) = (
    cbstate = (history = Read(), hyperparams = Read()),)

function on(::StepEnd, phase, logger::LogHyperParams, learner)
    history = learner.cbstate.history[phase]
    hyperparams = learner.cbstate.hyperparams
    for hparam in keys(hyperparams)
        val = last(last(hyperparams, hparam))
        log_to(
            logger.backends,
            Loggables.Value(val),
            string(hparam),
            history.steps,
            group = ("Step", "HParams"))
    end
end

# TODO: add support for logging histograms of layer activations and gradients
"""
    LogHistograms(backends...[; freq = 100]) <: Callback

Callback that logs histograms of model weights to [`LoggerBackend`](#)s
`backends` every `freq` steps.

If histograms should be logged every step, pass `freq = nothing`
"""
struct LogHistograms <: Callback
    backends::Tuple
    function LogHistograms(backends...; freq = 100)
        if isnothing(freq)
            return new(backends)
        else
            return throttle(new(backends), StepEnd, freq = freq)
        end
    end
end

stateaccess(::LogHistograms) = (model = Read(), cbstate = (history = Read(),))

function on(::StepEnd, phase::AbstractTrainingPhase, logger::LogHistograms, learner)
    history = learner.cbstate.history[phase]
    log_parameters(
        logger.backends,
        learner.model,
        "Model",
        history.steps,
        group = ("Step", "Histograms", "Weights"))
end

function log_parameters(backends, x, name, epochs; group)
    params = Flux.trainable(x)
    if isempty(params) && x isa AbstractArray
       log_to(
            backends,
            Loggables.Histogram(vec(x)),
            name,
            epochs,
            group = group)
    else
        for (pname, pval) in pairs(params)
            log_parameters(
                backends,
                pval,
                "$name.$pname",
                epochs,
                group = group)
        end
    end
end

"""
    LogVisualization(visfn, backends...[; freq = 100])

Logs images created by `visfn(learner.step)` to `backends` every `freq` steps.
"""
struct LogVisualization <: Callback
    visfn
    backends::Tuple
    function LogVisualization(visfn, backends...; freq = 100)
        cb = new(visfn, backends)
        if isnothing(freq)
            return cb
        else
            return throttle(cb, StepEnd, freq = freq)
        end
    end

end

stateaccess(::LogVisualization) = (step = Read(), cbstate = (history = Read(),))

function on(::StepEnd, phase::AbstractTrainingPhase, logger::LogVisualization, learner)
    history = learner.cbstate.history[phase]
    image = logger.visfn(learner.step)

    log_to(
        logger.backends,
        Loggables.Image(image),
        "Visualization",
        history.steps,
        group = ("Step",))
end


text:
using FileIO, JLD2
using Statistics, StatsBase
using PanGraph
using CairoMakie, ColorSchemes

delta(X,x,len) = minimum(mod.(abs.(X .- x), len))

function compare(qry, ref; cutoff=150)
    data = Int[]
    len  = length(Graphs.sequence(qry))

    # sample empirical distribution
    for i in 1:length(qry.node)
        x, y  = qry.position[i], qry.position[i+1]
        depth = Graphs.depth(qry.node[i].block)
        if depth ≤ cutoff && length(qry.node[i]) ≥ 1000
            push!(data, delta(ref.position, x, len))
            push!(data, delta(ref.position, y, len))
        end
    end

    # generate null distribution (uniform sampling of fixed number of breakpoints)
    x = sort(sample(1:len, length(data)÷2, replace=false))
    null = Array{Int,1}(undef,length(x))

    for i in 1:length(x)
        null[i] = delta(ref.position, x[i], len)
    end

    return (data, null)
end

function compare(graph)
    data = Int[]
    null = Int[]

    for (name, qry) in graph.pang.sequence
        d, n = compare(qry, graph.panx.sequence[name])
        append!(data,d); append!(null,n)
    end

    return data, null
end

function unmarshal(path)
    return (
        pang = open(Graphs.unmarshal,path.pang),
        panx = open(Graphs.unmarshal,path.panx),
    )
end

function group(args)
    i = 1
    sorted = sort(args)
    species = Dict(
        "kleb"  => "Klebsiella Pneumoniae",
        "myco"  => "Mycobacterium Tuberculosis",
        "proc"  => "Prochlorococcus Mariunus",
        "ecoli" => "Escherichia Coli",
        "helio" => "Helicobacter Pylori",
    )

    pairs = NamedTuple{(:pang,:panx,:name),Tuple{String,String,String}}[]
    while i ≤ length(sorted)
        push!(pairs, (
            pang=sorted[i],
            panx=sorted[i+1],
            name=species[basename(dirname(sorted[i]))],
        ))
        i += 2
    end

    return pairs
end

function shorten(name)
    words = split(name)

    return "$(uppercase(words[1][1])). $(words[2])"
end

function main(paths)
    fig = Figure(font="Latin Modern Math", fontsize=26)
    axis = Axis(fig[1,1],
        xlabel="Species",
        ylabel="distance to nearest gene (bp)",
        xticks=(1:5, map((p)->shorten(p.name), paths)),
        xticklabelrotation=π/6,
        yticks=(0:4, [L"10^0", L"10^1", L"10^2", L"10^3", L"10^4"]),
    )
    colors = cgrad(:Spectral_5, length(paths), categorical=true)

    for (i,path) in enumerate(paths)
        plots!(axis, path, i, colors[i])
    end

    fig
end

function plots!(axis, path, i, color)
    exported = "$(dirname(path.pang))/export.jld2"
    data, null = if !isfile(exported)
        data, null = path |> unmarshal |> compare
        FileIO.save(exported, Dict("data"=>data, "null"=>null))
        data, null
    else
        load(exported, "data", "null")
    end

    x = vcat(fill(i,length(data)+length(null)))
    y = log10.(vcat(data, null) .+ 1)
    s = vcat(fill(:left, size(data)), fill(:right, size(null)))
    c = vcat(fill((color,1.0),size(data)), fill((color,0.5),size(null)))

    violin!(axis, x, y, side=s, color=c)
    # cdfplot!(data.+1;
    #     color=color,
    #     linewidth=2,
    #     label="$(path.name) data"
    # )
    # cdfplot!(null.+1;
    #     color=color,
    #     linewidth=2,
    #     linestyle=:dashdot,
    #     label="$(path.name) null",
    # )
end

function save(plt)
    CairoMakie.save("figs/panx-compare.png", plt, px_per_unit=2)
end

if abspath(PROGRAM_FILE) == @__FILE__
    ARGS |> group |> main |> save
end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
@testset "Update boundaries with vertices" begin for b2index=1:16 begin for b1index = 1:16 for evpos = 1:4 begin begin
    if b1index == b2index
      continue
    end

    vpos = QuadTreeMeshes.POS(evpos)

    #print_with_color(:yellow, "------------------Testing ($(b1index), $(b2index), $vpos)\n")
    # add boundaries to single element
    x0, y0 = 2.0, 3.0
    dx, dy = 4.0, 4.0
    bb = GeometryTypes.SimpleRectangle(x0, y0, dx, dy)
    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)

    # add center vertex (center of child bounding box)
    elBB = QuadTreeMeshes.get_element_bounding_box(qt, 1)
    vx = elBB.x + 0.25 * elBB.w + ((vpos == QuadTreeMeshes.northEast || vpos == QuadTreeMeshes.southEast) ? 1 : 0) * 0.5 * elBB.w
    vy = elBB.y + 0.25 * elBB.h + ((vpos == QuadTreeMeshes.northWest || vpos == QuadTreeMeshes.northEast) ? 1 : 0) * 0.5 * elBB.h
    vPoint = QuadTreeMeshes.Point(vx, vy)

    push!(qt.vertices, vPoint)
    vindex = length(qt.vertices)

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

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

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

    # compute intersection points of boundary with bounding box
    b1, b2 = qt.vertices[vb1Index], qt.vertices[vb2Index]
    v = qt.vertices[vindex]
    b1vs = GeometryTypes.LineSegment(b1, v)
    vb2s = GeometryTypes.LineSegment(v, b2)

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

    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
end
end

@testset "Triangulate boundary template north south" begin
    bb = GeometryTypes.SimpleRectangle(2.0, 3.0, 4.0, 4.0)
    mesh = QuadTreeMeshes.QuadTreeMesh(bb)
    q = mesh.quadtree

    # add north-south vertices
    push!(q.vertices, QuadTreeMeshes.Point(2.1, 7.0))
    vnIndex = length(q.vertices)
    push!(q.vertices, QuadTreeMeshes.Point(2.7, 3.0))
    vsIndex = length(q.vertices)

    # triangulate NW element
    QuadTreeMeshes.triangulate_boundary_leave(mesh, 1, vnIndex, vsIndex)
    @test size(mesh.triangles) == (4,)

    # get triangulation
    @test !isnull(q.values[1])
    mesh_element = get(q.values[1])
    @test mesh_element.triangle_indices == [1,2,3,4]
    @test mesh.triangles[1].vertex_indices == FixedSizeArrays.Vec{3, QuadTreeMeshes.vertex_index}([1,6,3])
    @test mesh.triangles[2].vertex_indices == FixedSizeArrays.Vec{3, QuadTreeMeshes.vertex_index}([6,5,3])
    @test mesh.triangles[3].vertex_indices == FixedSizeArrays.Vec{3, QuadTreeMeshes.vertex_index}([6,2,5])
    @test mesh.triangles[4].vertex_indices == FixedSizeArrays.Vec{3, QuadTreeMeshes.vertex_index}([5,2,4])

    # plot triangulation
    filename = "triangulate_bndy_leave_ns.svg"
    Plots.plot(mesh)
    Plots.savefig(filename)
end

@testset "Triangulate template 0100" begin
    bb = GeometryTypes.SimpleRectangle(2.0, 3.0, 4.0, 4.0)
    mesh = QuadTreeMeshes.QuadTreeMesh(bb)
    q = mesh.quadtree

    # subdivide NE twice
    QuadTreeMeshes.subdivide!(q, 1, subdividefunc)
    firstElement = get_neighbour_from_index(q, 1, 2)
    QuadTreeMeshes.subdivide!(q, get(firstElement), subdividefunc)

    # triangulate NW element
    ne_element = get(get_neighbour_from_index(q, 1, 1))
    QuadTreeMeshes.triangulate_leave(mesh, ne_element)
    @test size(mesh.triangles) == (3,)

    # get triangulation
    @test !isnull(q.values[ne_element])
    mesh_element = get(q.values[ne_element])
    @test mesh_element.triangle_indices == [1,2,3]
    @test mesh.triangles[1].vertex_indices == FixedSizeArrays.Vec{3, QuadTreeMeshes.vertex_index}([8,5,13])
    @test mesh.triangles[2].vertex_indices == FixedSizeArrays.Vec{3, QuadTreeMeshes.vertex_index}([8,13,3])
    @test mesh.triangles[3].vertex_indices == FixedSizeArrays.Vec{3, QuadTreeMeshes.vertex_index}([3,13,6])

    # plot triangulation
    #filename = "triangulate_leave_0100.svg"
    #Plots.plot(mesh)
    #Plots.savefig(filename)
end

@testset "Triangulate template 1000" begin
    bb = GeometryTypes.SimpleRectangle(2.0, 3.0, 4.0, 4.0)
    mesh = QuadTreeMeshes.QuadTreeMesh(bb)
    q = mesh.quadtree

    # subdivide NW twice
    QuadTreeMeshes.subdivide!(q, 1, subdividefunc)
    firstElement = get_neighbour_from_index(q, 1, 1)
    QuadTreeMeshes.subdivide!(q, get(firstElement), subdividefunc)

    # triangulate NE element
    ne_element = get(get_neighbour_from_index(q, 1, 2))
    QuadTreeMeshes.triangulate_leave(mesh, ne_element)
    @test size(mesh.triangles) == (3,)

    # get triangulation
    @test !isnull(q.values[ne_element])
    mesh_element = get(q.values[ne_element])
    @test mesh_element.triangle_indices == [1,2,3]
    @test mesh.triangles[1].vertex_indices == FixedSizeArrays.Vec{3, QuadTreeMeshes.vertex_index}([5,9,14])
    @test mesh.triangles[2].vertex_indices == FixedSizeArrays.Vec{3, QuadTreeMeshes.vertex_index}([9,14,4])
    @test mesh.triangles[3].vertex_indices == FixedSizeArrays.Vec{3, QuadTreeMeshes.vertex_index}([14,4,6])

    # plot triangulation
    #filename = "triangulate_leave_1000.svg"
    #Plots.plot(mesh)
    #Plots.savefig(filename)
end

@testset "Triangulate template 0010" begin
    bb = GeometryTypes.SimpleRectangle(2.0, 3.0, 4.0, 4.0)
    mesh = QuadTreeMeshes.QuadTreeMesh(bb)
    q = mesh.quadtree

    # subdivide NE twice
    QuadTreeMeshes.subdivide!(q, 1, subdividefunc)
    firstElement = get_neighbour_from_index(q, 1, 2)
    QuadTreeMeshes.subdivide!(q, get(firstElement), subdividefunc)

    # triangulate SE element
    ne_element = get(get_neighbour_from_index(q, 1, 4))
    QuadTreeMeshes.triangulate_leave(mesh, ne_element)
    @test size(mesh.triangles) == (3,)

    # get triangulation
    @test !isnull(q.values[ne_element])
    mesh_element = get(q.values[ne_element])
    @test mesh_element.triangle_indices == [1,2,3]
    @test mesh.triangles[1].vertex_indices == FixedSizeArrays.Vec{3, QuadTreeMeshes.vertex_index}([7,12,5])
    @test mesh.triangles[2].vertex_indices == FixedSizeArrays.Vec{3, QuadTreeMeshes.vertex_index}([7,12,2])
    @test mesh.triangles[3].vertex_indices == FixedSizeArrays.Vec{3, QuadTreeMeshes.vertex_index}([2,12,9])

    # plot triangulation
    #filename = "triangulate_leave_0010.svg"
    #Plots.plot(mesh)
    #Plots.savefig(filename)
  end

  @testset "Triangulate template 0001" begin
    bb = GeometryTypes.SimpleRectangle(2.0, 3.0, 4.0, 4.0)
    mesh = QuadTreeMeshes.QuadTreeMesh(bb)
    q = mesh.quadtree

    # subdivide SE twice
    QuadTreeMeshes.subdivide!(q, 1, subdividefunc)
    firstElement = get_neighbour_from_index(q, 1, 4)
    QuadTreeMeshes.subdivide!(q, get(firstElement), subdividefunc)

    # triangulate NE element
    ne_element = get(get_neighbour_from_index(q, 1, 2))
    QuadTreeMeshes.triangulate_leave(mesh, ne_element)
    @test size(mesh.triangles) == (3,)

    # get triangulation
    @test !isnull(q.values[ne_element])
    mesh_element = get(q.values[ne_element])
    @test mesh_element.triangle_indices == [1,2,3]
    @test mesh.triangles[1].vertex_indices == FixedSizeArrays.Vec{3, QuadTreeMeshes.vertex_index}([5,11,6])
    @test mesh.triangles[2].vertex_indices == FixedSizeArrays.Vec{3, QuadTreeMeshes.vertex_index}([6,11,4])
    @test mesh.triangles[3].vertex_indices == FixedSizeArrays.Vec{3, QuadTreeMeshes.vertex_index}([11,9,4])

    # plot triangulation
    #filename = "triangulate_leave_0001.svg"
    #Plots.plot(mesh)
    #Plots.savefig(filename)
  end

  @testset "Triangulate template 0000" begin
    bb = GeometryTypes.SimpleRectangle(2.0, 3.0, 4.0, 4.0)
    mesh = QuadTreeMeshes.QuadTreeMesh(bb)
    q = mesh.quadtree

    # triangulate first and only element
    QuadTreeMeshes.triangulate_leave(mesh, 1)
    @test size(mesh.triangles) == (2,)

    # get triangulation
    @test !isnull(q.values[1])
    mesh_element = get(q.values[1])
    @test mesh_element.triangle_indices == [1,2]
    @test mesh.triangles[1].vertex_indices == FixedSizeArrays.Vec{3, QuadTreeMeshes.vertex_index}([1,2,3])
    @test mesh.triangles[2].vertex_indices == FixedSizeArrays.Vec{3, QuadTreeMeshes.vertex_index}([2,4,3])

    # plot triangulation
    #filename = "triangulate_leave_0000.svg"
    #Plots.plot(mesh)
    #Plots.savefig(filename)
  end


text:
@testset "Insertion" begin
    t = SegmentMap{Int,Int}()
    n = 100000

    @testset "random insertions" begin
        for v in 1:n
            k = randsegment()
            t[k] = v
        end
        @test issorted(collect(keys(t)))
        @test validkeys(t)
        @test validparents(t)
        @test validsiblings(t)
    end

    @testset "ordered insertions" begin
        for v in 1:n
            k = (v, rand(v:maxend))
            t[k] = v
        end
        @test issorted(collect(keys(t)))
        @test validkeys(t)
        @test validparents(t)
        @test validsiblings(t)
    end

    @testset "reverse ordered insertions" begin
        for v in 1:n
            k = (maxend - v, rand((maxend - v):maxend))
            t[k] = v
        end
        @test issorted(collect(keys(t)))
        @test validkeys(t)
        @test validparents(t)
        @test validsiblings(t)
    end

    @testset "bulk insertion" begin
        segments = SegmentValue{Int,Int}[]
        for v in 1:n
            a, b = randsegment()
            push!(segments, SegmentValue{Int,Int}(a, b, v))
        end
        sort!(segments)
        t = SegmentMap{Int,Int}(segments)

        @test collect(keys(t)) == [Segment{Int}(segment.first, segment.last)
                                   for segment in segments]
        @test collect(values(t)) == [segment.value for segment in segments]
        @test issorted(collect(keys(t)))
        @test validkeys(t)
        @test validparents(t)
        @test validsiblings(t)

        # don't break on empty arrays
        t = SegmentMap{Int,Int}(SegmentValue{Int,Int}[])
        @test issorted(collect(keys(t)))
        @test validkeys(t)
        @test validparents(t)
        @test validsiblings(t)

        @test_throws MethodError SegmentTree{Int,Int}(segments, Int[1])
        shuffle!(segments)
        @test_throws UndefVarError SegmentTree{Int,Int}(segments, vals)
    end
end

@testset "Updates" begin
    t = SegmentMap{Int,Int}()
    n = 100000
    ks = [randsegment() for _ in 1:n]
    for (v, k) in enumerate(ks)
        t[k] = v
    end

    shuffle!(ks)
    for (v, k) in enumerate(ks)
        t[k] = v
    end

    @test issorted(collect(keys(t)))
    @test validkeys(t)
    @test validparents(t)
    @test validsiblings(t)
end

@testset "Deletion" begin
    # Set B to smaller number to try to induce more weird merges and borrows and
    # such
    t = SegmentBTree{Int,SegmentValue{Int,Int},8}()
    n = 100000

    # insert n random segments
    segments = [randsegment() for _ in 1:n]
    for (i, segment) in enumerate(segments)
        t[segment] = i
    end
    shuffle!(segments)

    # delete non-existant is a nop
    @test delete!(t, (-1, -1)).n === n

    # delete one
    segment = pop!(segments)
    delete!(t, segment)
    @test validkeys(t)
    @test validparents(t)
    @test validsiblings(t)
    @test issorted(collect(keys(t)))
    @test !haskey(t, segment)
    @test all(segment -> haskey(t, segment), segments)

    # delete a random 50%
    for segment in 1:(n - 1) / 2
        segment = pop!(segments)
        delete!(t, segment)
    end
    @test validkeys(t)
    @test validparents(t)
    @test validsiblings(t)
    @test issorted(collect(keys(t)))
    @test all(segment -> haskey(t, segment), segments)

    # delete the rest
    for segment in segments
        delete!(t, segment)
    end
    @test isempty(t)
    @test validkeys(t)
    @test validparents(t)
    @test validsiblings(t)

    # delete left-to-right
    shuffle!(segments)
    for (i, segment) in enumerate(segments)
        t[segment] = i
    end
    sort!(segments)
    for segment in segments
        delete!(t, segment)
    end
    @test isempty(t)
    @test validkeys(t)
    @test validparents(t)
    @test validsiblings(t)

    # delete right-to-left
    shuffle!(segments)
    for (i, segment) in enumerate(segments)
        t[segment] = i
    end
    reverse!(segments)
    for segment in segments
        delete!(t, segment)
    end
    @test isempty(t)
    @test validkeys(t)
    @test validparents(t)
    @test validsiblings(t)

    # delete from the middle
    shuffle!(segments)
    for (i, segment) in enumerate(segments)
        t[segment] = i
    end
    mid = div(length(segments), 2)
    for i in mid + 1:length(segments)
        delete!(t, segments[i])
    end
    for i in mid:-1:1
    delete!(t, segments[i])
    end
    @test isempty(t)
    @test validkeys(t)
    @test validparents(t)
    @test validsiblings(t)

    # contrive a special case: merge with the left when we have a non-null
    # right
    t = SegmentBTree{Int,SegmentValue{Int,Int},4}()
    for i in 1:24
        t[(i, i)] = i
    end
    @test isa(t.root, InternalNode)
    @test length(t.root.children) === 3
    @test length(t.root.children[2].children) === 2
    delete!(t, (15, 15))
    keys_to_delete = collect(keys(t))
    for key in keys_to_delete
        delete!(t, key)
    end
    @test isempty(t)
    @test validkeys(t)
    @test validparents(t)
    @test validsiblings(t)
end


text:
"""
QQQQ
"""
function maph_init_test()
    # desired_absorbptions = [0.3,0.5,0.1,0.08,0.02]
    # desired_means = [2.0,3.1,2.3,4.5,0.2]
    # out = MAPHDist(20,desired_absorbptions,desired_means,[1.1,2.3,2.0,0.33,1.5])
    # # @show model_size(out)
    # # data = [rand(out) for _ in 1:10^5]
    # # pi_est = [count((x)->x.a == i, data)/length(data) for i in 1:5]
    # # μ_est = [mean(first.(filter((x)->x.a == i, data))) for i in 1:5]
    # (absorption_probs(out) ≈ desired_absorbptions') || return false

    Λ₄, λ45, λ54, Λ₅ = 20, 17, 7, 10
    μ41, μ42, μ43, μ51, μ52, μ53 = 1, 1, 1, 1, 1, 1
    T_example = [-Λ₄ λ45; λ54 -Λ₅]
    T0_example = [μ41 μ42 μ43; μ51 μ52 μ53]
    initial_dist = [0.5,0.5]
    maph = MAPHDist(initial_dist', 3*T_example, 3*T0_example)
    @show absorption_probs(maph)
    @show mean(maph)

    # mean(out) ≈ desired_means || return false
    # @show pi_est
    # @show μ_est
    # @show pi_est .* μ_est
    return true
end


text:
#
# Optimizing pipelines.
#

struct RewriteMemo
    op_map::Dict{Any,Dict{Vector{Any},Pipeline}}

    RewriteMemo() = new(Dict{Any,Dict{Vector{Any},Pipeline}}())
end

function (memo::RewriteMemo)(p::Pipeline)
    args_map = get!(memo.op_map, p.op) do
        Dict{Vector{Any},Pipeline}()
    end
    get!(args_map, p.args) do
        p
    end
end

function (memo::RewriteMemo)(ps::Vector{Pipeline})
    if isempty(ps)
        memo(pass())
    elseif length(ps) == 1
        memo(ps[1])
    else
        memo(chain_of(ps))
    end
end

function unchain(p)
    if p.op === pass
        return Pipeline[]
    elseif p.op === chain_of
        return collect(Pipeline, p.args[1])
    else
        return Pipeline[p]
    end
end

function rewrite_all(p::Pipeline; memo=RewriteMemo())::Pipeline
    rewrite_common(rewrite_simplify(p, memo=memo), memo=memo)
end

@inline function rewrite_with(f, memo, p)
    @match_pipeline if (p ~ chain_of(qs))
        chain = Pipeline[]
        for q in qs
            append!(chain, unchain(f(memo, q)))
        end
        return memo(chain)
    end
    args = Any[f(memo, arg) for arg in p.args]
    memo(Pipeline(p.op, args=args))
end

#
# Local simplification.
#

function rewrite_simplify(p::Pipeline; memo=RewriteMemo())::Pipeline
    rewrite_simplify(memo, p) |> designate(p.sig)
end

function rewrite_simplify(memo::RewriteMemo, p::Pipeline)
    @match_pipeline if (p ~ chain_of(_))
        chain = Pipeline[]
        simplify_and_push!(memo, chain, p)
        return memo(chain)
    end
    p′ = rewrite_with(rewrite_simplify, memo, p)
    # with_column(N, pass()) => pass()
    @match_pipeline if (p′ ~ with_column(_, pass()))
        p′ =  memo(pass())
    # with_elements(pass()) => pass()
    elseif (p′ ~ with_elements(pass()))
        p′ = memo(pass())
    end
    p′
end

rewrite_simplify(memo::RewriteMemo, ps::Vector{Pipeline}) =
    Pipeline[rewrite_simplify(memo, p) for p in ps]

rewrite_simplify(memo::RewriteMemo, @nospecialize other) = other

function simplify_and_push!(memo::RewriteMemo, chain::Vector{Pipeline}, p::Pipeline)
    top = !isempty(chain) ? chain[end] : memo(pass())
    @match_pipeline if (p ~ pass())
    elseif (p ~ chain_of(qs))
        for q in qs
            if (q ~ chain_of(_))
                simplify_and_push!(memo, chain, q)
            else
                simplify_and_push!(memo, chain, rewrite_simplify(memo, q))
            end
        end
    elseif (p ~ tuple_of(lbls, cols))
        head = nothing
        for col in cols
            if (col ~ chain_of([q, _...]))
                col = q
            end
            if (col ~ filler(_)) || (col ~ block_filler(_, _)) || (col ~ null_filler())
            elseif (col ~ column(k)) && (head === nothing || head === col)
                head = col
            else
                head = missing
                break
            end
        end
        if head === missing
            push!(chain, p)
        elseif head === nothing
            empty!(chain)
            push!(chain, p)
        else
            cols′ = Pipeline[]
            for col in cols
                if col === head
                    col = memo(pass())
                elseif (col ~ chain_of([q, qs...]))
                    if q === head
                        col = memo(qs)
                    end
                end
                push!(cols′, col)
            end
            simplify_and_push!(memo, chain, head)
            simplify_and_push!(memo, chain, memo(tuple_of(lbls, cols′)))
        end
    # chain_of(p, filler(val)) => filler(val)
    elseif (p ~ filler(_)) || (p ~ block_filler(_, _)) || (p ~ null_filler())
        empty!(chain)
        push!(chain, p)
    # chain_of(tuple_of(p1, ..., pn), column(k)) => pk
    elseif (top ~ tuple_of(lbls, cols)) && (p ~ column(lbl))
        k = lbl isa Symbol ? findfirst(isequal(lbl), lbls) : lbl
        qs = unchain(cols[k])
        pop!(chain)
        for q in qs
            simplify_and_push!(memo, chain, q)
        end
    # chain_of(tuple_of(..., pk, ...), with_column(k, q)) => tuple_of(..., chain_of(pk, q), ...)
    elseif (top ~ tuple_of(lbls, cols)) && (p ~ with_column(lbl, q))
        k = lbl isa Symbol ? findfirst(isequal(lbl), lbls) : lbl
        pop!(chain)
        qs = unchain(cols[k])
        simplify_and_push!(memo, qs, q)
        cols′ = copy(cols)
        cols′[k] = memo(qs)
        simplify_and_push!(memo, chain, memo(tuple_of(lbls, cols′)))
    # chain_of(tuple_of(chain_of(p, wrap()), ...), tuple_lift(f)) => chain_of(tuple_of(p, ...), tuple_lift(f))
    elseif (top ~ tuple_of(lbls, cols)) && (p ~ tuple_lift(_))
        cols′ = Pipeline[]
        for col in cols
            qs = unchain(col)
            while (qs ~ [_..., wrap()])
                pop!(qs)
            end
            push!(cols′, memo(qs))
        end
        pop!(chain)
        push!(chain, memo(tuple_of(lbls, cols′)))
        push!(chain, p)
    # chain_of(with_column(k, wrap()), distribute(k)) => wrap()
    elseif (top ~ with_column(j, wrap())) && (p ~ distribute(k)) && j == k
        pop!(chain)
        simplify_and_push!(memo, chain, memo(wrap()))
    # chain_of(with_column(k, chain_of(p, wrap())), distribute(k)) => chain_of(with_column(k, p), wrap())
    elseif (top ~ with_column(j, chain_of([qs..., wrap()]))) && (p ~ distribute(k)) && j == k
        pop!(chain)
        push!(chain, memo(with_column(k, memo(qs))))
        push!(chain, memo(wrap()))
    # chain_of(wrap(), flatten()) => pass()
    elseif (top ~ wrap()) && (p ~ flatten())
        pop!(chain)
    # chain_of(wrap(), with_elements(p)) => chain_of(p, wrap())
    elseif (top ~ wrap()) && (p ~ with_elements(q))
        pop!(chain)
        qs = unchain(q)
        for q in qs
            simplify_and_push!(memo, chain, q)
        end
        simplify_and_push!(memo, chain, top)
    # chain_of(wrap(), lift(f)) => lift(f)
    elseif (top ~ wrap()) && (p ~ lift(_))
        pop!(chain)
        simplify_and_push!(memo, chain, p)
    # chain_of(with_elements(p), with_elements(q)) => with_elements(chain_of(p, q))
    elseif (top ~ with_elements(q1)) && (p ~ with_elements(q2))
        pop!(chain)
        qs = unchain(q1)
        for q in unchain(q2)
            simplify_and_push!(memo, qs, q)
        end
        simplify_and_push!(memo, chain, memo(with_elements(memo(qs))))
    # chain_of(with_elements(wrap()), flatten()) => pass()
    elseif (top ~ with_elements(wrap())) && (p ~ flatten())
        pop!(chain)
    # chain_of(with_elements(chain_of(p, wrap())), flatten()) => with_elements(p)
    elseif (top ~ with_elements(chain_of([qs..., wrap()]))) && (p ~ flatten())
        pop!(chain)
        simplify_and_push!(memo, chain, memo(with_elements(memo(qs))))
    # chain_of(flatten(), with_elements(column(k))) => chain_of(with_elements(with_elements(column(k))), flatten())
    elseif (top ~ flatten()) && (p ~ with_elements(column(k)))
        pop!(chain)
        simplify_and_push!(memo, chain, memo(with_elements(p)))
        simplify_and_push!(memo, chain, top)
   # chain_of(flatten(), with_elements(chain_of(column(k), p))) => chain_of(with_elements(with_elements(column(k))), flatten(), with_elements(p))
    elseif (top ~ flatten()) && (p ~ with_elements(chain_of([column(k), qs...])))
        pop!(chain)
        simplify_and_push!(memo, chain, memo(with_elements(memo(with_elements(memo(column(k)))))))
        simplify_and_push!(memo, chain, top)
        simplify_and_push!(memo, chain, memo(with_elements(memo(qs))))
    # chain_of(distribute(k), with_elements(column(k))) => column(k)
    elseif (top ~ distribute(j)) && (p ~ with_elements(column(k))) && j == k
        pop!(chain)
        simplify_and_push!(memo, chain, memo(column(k)))
    # chain_of(distribute(k), with_elements(chain_of(column(k), p))) => chain_of(column(k), with_elements(p))
    elseif (top ~ distribute(j)) && (p ~ with_elements(chain_of([column(k), qs...]))) && j == k
        pop!(chain)
        simplify_and_push!(memo, chain, memo(column(k)))
        simplify_and_push!(memo, chain, memo(with_elements(memo(qs))))
    # chain_of(sieve_by(), with_elements(column(k))) => chain_of(with_column(1, column(k)), sieve_by())
    elseif (top ~ sieve_by()) && (p ~ with_elements(column(k)))
        pop!(chain)
        simplify_and_push!(memo, chain, memo(with_column(1, memo(column(k)))))
        simplify_and_push!(memo, chain, top)
    else
        push!(chain, p)
    end
    nothing
end

#
# Common subexpression elimination.
#

function rewrite_common(p::Pipeline; memo=RewriteMemo())::Pipeline
    rewrite_common(memo, p) |> designate(p.sig)
end

function rewrite_common(memo::RewriteMemo, p::Pipeline)
    @match_pipeline if (p ~ tuple_of(_, _))
        chain = Pipeline[]
        l, r = pull_common(memo, p)
        append!(chain, unchain(rewrite_common(memo, l)))
        push!(chain, rewrite_with(rewrite_common, memo, r))
        memo(chain)
    else
        rewrite_with(rewrite_common, memo, p)
    end
end

rewrite_common(memo::RewriteMemo, ps::Vector{Pipeline}) =
    Pipeline[rewrite_common(memo, p) for p in ps]

rewrite_common(memo::RewriteMemo, @nospecialize other) = other

function pull_common(memo::RewriteMemo, p::Pipeline)
    l = memo(pass())
    r = p
    @match_pipeline if (p ~ tuple_of(lbls, cols)) && length(cols) > 1
        i = l
        deps = Dict{Pipeline,Vector{Pipeline}}(i => Pipeline[], p => cols)
        seen = Set{Pipeline}()
        dups = Set{Pipeline}()
        for col in cols
            parts = Pipeline[i]
            collect_parts!(memo, parts, deps, i, col)
            for q in parts
                if q in seen
                    push!(dups, q)
                end
            end
            for q in parts
                push!(seen, q)
            end
        end
        if length(dups) > 1
            seen = Set{Pipeline}()
            basis = Pipeline[]
            stack = Pipeline[p]
            while !isempty(stack)
                q = pop!(stack)
                for dep in deps[q]
                    if !(dep in seen)
                        if dep in dups
                            push!(basis, dep)
                        else
                            push!(stack, dep)
                        end
                        push!(seen, dep)
                    end
                end
            end
            w = length(basis)
            if w > 1
                basis = Pipeline[q for q in basis if !(q ~ column(_))]
                if length(basis) < w && !(i in basis)
                    push!(basis, i)
                end
            end
            if !(basis ~ [pass()])
                repl = Dict{Pipeline,Pipeline}()
                if length(basis) == 1
                    l = basis[1]
                    repl[basis[1]] = i
                else
                    l_cols = Pipeline[]
                    for (k, b) in enumerate(basis)
                        push!(l_cols, b)
                        repl[b] = memo(column(k))
                    end
                    l = tuple_of(Symbol[], l_cols)
                end
                r = replace_parts(memo, repl, i, p)
                @assert r !== nothing p
            end
        end
    end
    (l, r)
end

function collect_parts!(memo, parts, deps, i, p)
    @match_pipeline if (p ~ pass())
        o = i
    elseif (p ~ chain_of(qs))
        o = i
        for q in qs
            o = collect_parts!(memo, parts, deps, o, q)
        end
    elseif (p ~ tuple_of(lbls, cols))
        os = Pipeline[]
        for col in cols
            push!(os, collect_parts!(memo, parts, deps, i, col))
        end
        chain = unchain(i)
        push!(chain, p)
        o = memo(chain)
        push!(parts, o)
        get!(deps, o, os)
    elseif (p ~ with_elements(chain_of(qs)))
        chain = unchain(i)
        qs′ = Pipeline[]
        for q in qs
            push!(qs′, q)
            chain′ = copy(chain)
            push!(chain′, memo(with_elements(memo(copy(qs′)))))
            o = memo(chain′)
            get!(deps, o, Pipeline[i])
            i = o
        end
    else
        chain = unchain(i)
        push!(chain, p)
        o = memo(chain)
        push!(parts, o)
        if (p ~ filler(_)) || (p ~ block_filler(_, _)) || (p ~ null_filler())
            get!(deps, o, Pipeline[])
        else
            get!(deps, o, Pipeline[i])
        end
    end
    o
end

function replace_parts(memo, repl, i, p)
    o = memo(vcat(unchain(i), unchain(p)))
    if o in keys(repl)
        return repl[o]
    end
    @match_pipeline if (p ~ chain_of(qs))
        chain′ = unchain(i)
        append!(chain′, qs)
        for k in reverse(eachindex(qs))
            pop!(chain′)
            i′ = memo(copy(chain′))
            o′ = replace_parts(memo, repl, i′, qs[k])
            if o′ !== nothing
                chain = unchain(o′)
                append!(chain, qs[k+1:end])
                return memo(chain)
            end
        end
    elseif (p ~ tuple_of(lbls, cols))
        cols′ = Pipeline[]
        for col in cols
            o = replace_parts(memo, repl, i, col)
            if o !== nothing
                push!(cols′, o)
            end
        end
        if length(cols′) == length(cols)
            return tuple_of(lbls, cols′)
        end
    elseif (p ~ with_elements(chain_of(qs)))
        for k in lastindex(qs)-1:-1:firstindex(qs)+1
            chain′ = unchain(i)
            push!(chain′, memo(with_elements(memo(qs[1:k-1]))))
            i′ = memo(chain′)
            if i′ in keys(repl)
                o′ = repl[i′]
                chain = unchain(o′)
                push!(chain, memo(with_elements(memo(qs[k+1:end]))))
                return memo(chain)
            end
        end
    elseif (p ~ filler(_)) || (p ~ block_filler(_, _)) || (p ~ null_filler())
        return p
    end
    if i in keys(repl)
        return memo(vcat(unchain(repl[i]), unchain(p)))
    end
    nothing
end

function rewrite_linearize(p::Pipeline; memo=RewriteMemo())::Pipeline
    chain_of(linearize(p)) |> designate(signature(p))
end

function linearize(p::Pipeline)::Vector{Pipeline}
    retval = Pipeline[]
    @match_pipeline if (p ~ pass())
        nothing
    elseif (p ~ chain_of(qs))
        for q in qs
            append!(retval, linearize(q))
        end
    elseif (p ~ with_elements(q))
        for r in linearize(q)
            push!(retval, with_elements(r))
        end
    elseif (p ~ with_column(lbl, q))
        for r in linearize(q)
            push!(retval, with_column(lbl, r))
        end
    elseif (p ~ tuple_of(lbls, cols::Vector{Pipeline}))
        push!(retval, tuple_of(lbls, length(cols)))
        for (idx, q) in enumerate(cols)
            for r in linearize(q)
                push!(retval, with_column(idx, r))
            end
        end
    else
        push!(retval, p)
    end
    return retval
end

function delinearize!(vp::Vector{Pipeline})::Pipeline
    chain = Vector{Pipeline}()
    while length(vp) > 0
        @match_pipeline if (vp[1] ~ with_elements(arg))
            push!(chain, delinearize_block!(vp))
        elseif (vp[1] ~ tuple_of(lbls, width::Int))
            push!(chain, delinearize_tuple!(vp, lbls, width))
        else
            push!(chain, popfirst!(vp))
        end
    end
    return chain_of(chain...)
end

function delinearize_block!(vp::Vector{Pipeline})::Pipeline
    chain = Vector{Pipeline}()
    @match_pipeline while (vp ~ [with_elements(p), _...])
        popfirst!(vp)
        push!(chain, p)
    end
    return with_elements(delinearize!(chain))
end

function delinearize_tuple!(vp::Vector{Pipeline}, lbls, width)::Pipeline
    popfirst!(vp) # drop the `tuple_of`
    slots = [Pipeline[] for x in 1:width]
    @match_pipeline while (vp ~ [with_column(idx, p), _...])
        popfirst!(vp)
        if !isa(idx, Int)
            idx = findfirst(==(idx), lbls)
            @assert !isnothing(idx)
        end
        push!(slots[idx], p)
    end
    return tuple_of(lbls, [delinearize!(cv) for cv in slots])
end


text:
module testMod1;

export test;

function test(x::Int64)
   sum([i^2 for i=1:x])
end

end

max_stars_repo_name: aniongithub/node-julia

max_issues_repo_name: aniongithub/node-julia

max_forks_repo_name: aniongithub/node-julia


text:
"""
	extrapolate(spl, xrange, extr::Symbol)

Extrapolate a spline to the limits given by `xrange`, using the extrapolation type `extr`.

- alias: `extrap`

Allowed modes for `extr`:
- `:zero`
- `:boundary`
- `:linear`
- `:quadratic`
- `:smooth` (smooth blending to zero for Hermite splines)
"""
function extrapolate(spl, xrange, extr::Symbol)
    envx = Float64.(collect(xrange))
    xs = 0.0 # start of extrapolatin range
    xe = 0.0 # end of extrapolation range
    envx[1] < spl.x[1] ? xs = envx[1] : xs = spl.x[1]
    envx[end] > spl.x[end] ? xe = envx[end] : xe = spl.x[end]

    x = copy(spl.x)
    a = copy(spl.a)
	nf = nextfloat(Float64(x[end]))

    if typeof(spl) in (LinearSpline, QuadraticSpline, CubicSpline, HermiteSpline)
        b = copy(spl.b)
    end
    if typeof(spl) in (QuadraticSpline, CubicSpline)
        c = copy(spl.c)
    end
    if typeof(spl) == CubicSpline
        d = copy(spl.d)
    end
	if extr == :zero
		if xs < x[1]
			if typeof(spl) == HermiteSpline
				pushfirst!(a, 0.0)
				pushfirst!(b, 0.0)
				pushfirst!(x, prevfloat(x[1]))
			end
			x[1] = prevfloat(x[1])
			pushfirst!(x, xs)
			pushfirst!(a, 0.0)
			if typeof(spl) in (LinearSpline, QuadraticSpline, CubicSpline, HermiteSpline)
				pushfirst!(b, 0.0)
			end
			if typeof(spl) in (QuadraticSpline,CubicSpline)
				pushfirst!(c, 0.0)
			end
			if typeof(spl) == CubicSpline
				pushfirst!(d, 0.0)
			end
		end
		if xe ≥ x[end]
		if typeof(spl) == HermiteSpline
			push!(x, nf)
			push!(x, xe)
			push!(a, 0.0)
			push!(a, 0.0)
			push!(b, 0.0)
			push!(b, 0.0)
		end
			if typeof(spl) in (LinearSpline, QuadraticSpline, CubicSpline)
				x[end] = nf
				a[end] = 0.0
				push!(x, xe)
				push!(a, 0.0)
				push!(b, 0.0)
			end
			if typeof(spl) in (QuadraticSpline, CubicSpline)
				push!(c, 0.0)
			end
			if typeof(spl) == CubicSpline
				push!(d, 0.0)
			end
		end
	elseif extr == :boundary
		de = x[end] - x[end-1]
		if xs < x[1]
			if typeof(spl) == HermiteSpline
				pushfirst!(a, a[1])
				pushfirst!(b, 0.0)
				pushfirst!(x, prevfloat(x[1]))
			end
			pushfirst!(x, xs)
			pushfirst!(a, a[1])
			if typeof(spl) in (LinearSpline, QuadraticSpline, CubicSpline, HermiteSpline)
				pushfirst!(b, 0.0)
			end
			if typeof(spl) in (QuadraticSpline, CubicSpline)
				pushfirst!(c, 0.0)
			end
			if typeof(spl) == CubicSpline
				pushfirst!(d, 0.0)
			end
		end
		if xe ≥ x[end]
			if typeof(spl) == LinearSpline
				push!(a, a[end])
				push!(b, 0.0)
			end
			if typeof(spl) == QuadraticSpline
				push!(a, a[end])
				push!(b, 0.0)
				push!(c, 0.0)
				x[end] = nf
			end
			if typeof(spl) == CubicSpline
				de2 = de * de
				push!(a,a[end])
				push!(b,0.0)
				push!(c,0.0)
				push!(d,0.0)
				x[end] = nf
			end
			if typeof(spl) == HermiteSpline
				push!(a, a[end])
				push!(a, a[end])
				push!(x, nf)
				push!(b, 0.0)
				push!(b, 0.0)
			end
		end
		push!(x,xe)
	elseif extr == :linear
		if xs < x[1]
			pushfirst!(x, xs)
			pushfirst!(a, a[1] + b[1] * (x[1] - x[2]))
			if typeof(spl) in (LinearSpline, QuadraticSpline, CubicSpline, HermiteSpline)
				pushfirst!(b, b[1])
			end
			if typeof(spl) in (QuadraticSpline,CubicSpline)
				pushfirst!(c, 0.0)
			end
			if typeof(spl) == CubicSpline
				pushfirst!(d, 0.0)
			end
		end
		if xe ≥ x[end]
			dn = xe - x[end]
			dn2 = dn * dn
			de = x[end] - x[end-1]
			de2 = de * de
			if typeof(spl) == LinearSpline
				push!(a, a[end] + b[end] * dn)
				push!(b, b[end])
			end
			if typeof(spl) == QuadraticSpline
				push!(b, b[end] + c[end] * dn)
				push!(a, a[end] + b[end] * dn)
				push!(c, 0.0)
			end
			if typeof(spl) == CubicSpline
				push!(b, 3 * d[end] * de2 + 2 * c[end] * de + b[end])
				push!(a, a[end] + b[end] * dn)
				push!(c, 0.0)
				push!(d, 0.0)
				x[end] = nf
			end
			if typeof(spl) == HermiteSpline
				push!(a, a[end] + b[end] * dn)
				push!(b, b[end])
			end
		end
		push!(x,xe)
    elseif extr == :quadratic
        if xs < x[1]
            pushfirst!(x, xs)
            if typeof(spl) in (LinearSpline, QuadraticSpline, CubicSpline)
                pushfirst!(b, x[1] * ((a[3] - a[2]) + x[2] * (a[1] - b[3]) + x[3] * (a[2] - a[1])) /
							((x[1] - x[2]) * (x[1] - x[3]) * (x[2] - x[3])))
                pushfirst!(a, a[1] / (b[1] * (x[1] - xs)))
            end
            if typeof(spl) in (QuadraticSpline, CubicSpline)
                pushfirst!(c, 0.0)
            end
            if typeof(spl) == CubicSpline
                pushfirst!(d, 0.0)
            end
        end
        if xe ≥ x[end]
            push!(x, xe)
            if typeof(spl) in (LinearSpline, QuadraticSpline, CubicSpline)
                push!(b, x[end-2] * ((a[end] - a[end-1]) + x[end-1] * (a[end-2] - b[end]) + x[end] * (a[end-1] - a[end-2])) /
						((x[end-2] - x[end-1]) * (x[end-2] - x[end]) * (x[end-1] - x[end])))
                push!(a, a[end]/(b[end]*(x[end]-xs)))
            end
            if typeof(spl) in (QuadraticSpline, CubicSpline)
                push!(c, 0.0)
            end
            if typeof(spl) == CubicSpline
                push!(d, 0.0)
            end
        end
	elseif extr == :smooth
		pushfirst!(x, xs)
		pushfirst!(x, prevfloat(xs))
		pushfirst!(a, 0.0)
		pushfirst!(a, 0.0)
		push!(x, xe)
		push!(x, nextfloat(xe))
		push!(a, 0.0)
		push!(a, 0.0)
		smooth = hermitespline(x, a, :kruger)
		a = smooth.a
		b = smooth.b
    else
		throw(DomainError(extr, "Extrapolation type does not exist."))
    end
    if typeof(spl) == NearestNeighborSpline
        NearestNeighborSpline(a, x)
    elseif typeof(spl) == LinearSpline
        LinearSpline(b, a, x)
    elseif typeof(spl) == QuadraticSpline
        QuadraticSpline(c, b, a, x)
	elseif typeof(spl) == CubicSpline
        CubicSpline(d, c, b, a, x)
	else
		HermiteSpline(b, a, x)
    end
end

"""
	extrapolate(spl, env = SpectralVis.SPECENV)

Extrapolate spline according to spectral environment.
"""
function extrapolate(spl, env = SpectralVis.SPECENV)
    extrapolate(spl, env.λmin:Δλ:λmax, env.ex)
end

# alias
extrap = extrapolate


text:
# Benchmark the counting of alanine residues in a PDB file

using BioStructures

pdb_filepath = "data/1AKE.pdb"
struc = read(pdb_filepath, PDB)

function counter()
    alanineselector(res::AbstractResidue) = resnameselector(res, ["ALA"])
    return countresidues(struc, alanineselector)
end

# Run to JIT compile
counter()

elapsed = @elapsed counter()

println(elapsed)

max_stars_repo_name: wojdyr/pdb-benchmarks

max_issues_repo_name: wojdyr/pdb-benchmarks

max_forks_repo_name: wojdyr/pdb-benchmarks


text:
"Summary of optimization using the NLopt package"
type OptSummary
    initial::Vector{Float64}
    final::Vector{Float64}
    fmin::Float64
    feval::Int
    geval::Int
    optimizer::Symbol
end
function OptSummary(initial::Vector{Float64},optimizer::Symbol)
    OptSummary(initial,initial,Inf,-1,-1,optimizer)
end

type LMM <: StatsBase.RegressionModel
    trms::Vector
    Λ::Vector
    A::Matrix        # symmetric cross-product blocks (upper triangle)
    R::Matrix        # right Cholesky factor in blocks.
    opt::OptSummary
end

function LMM(Rem::Vector, Λ::Vector, X::AbstractMatrix, y::Vector)
    all(x->isa(x,AbstractReMat),Rem) ||
        throw(ArgumentError("Elements of Rem should be AbstractReMat's"))
    all(x->isa(x,ParamLowerTriangular),Λ) ||
        throw(ArgumentError("Elements of Λ should be ParamLowerTriangular"))
    n,p = size(X)
    all(t -> size(t,1) == n,Rem) && length(y) == n || throw(DimensionMismatch("n not consistent"))
    nreterms = length(Rem)
    length(Λ) == nreterms && all(i->chksz(Rem[i],Λ[i]),1:nreterms) ||
        throw(DimensionMismatch("Rem and Λ"))
    nt = nreterms + 1
    trms = Array(Any,nt)
    for i in eachindex(Rem) trms[i] = Rem[i] end
    trms[end] = hcat(X,y)
    A = fill!(Array(Any,(nt,nt)),nothing)
    R = fill!(Array(Any,(nt,nt)),nothing)
    for j in 1:nt
        for i in 1:j
            A[i,j] = densify(trms[i]'trms[j])
            R[i,j] = copy(A[i,j])
        end
    end
    for j in 2:nreterms
        if isa(R[j,j],Diagonal) || isa(R[j,j],HBlkDiag)
            for i in 1:(j-1)     # check for fill-in
                if !isdiag(A[i,j]'A[i,j])
                    for k in j:nt
                        R[j,k] = full(R[j,k])
                    end
                end
            end
        end
    end
    LMM(trms,Λ,A,R,OptSummary(mapreduce(x->x[:θ],vcat,Λ),:None))
end

LMM(re::Vector,Λ::Vector,X::AbstractMatrix,y::DataVector) = LMM(re,Λ,X,convert(Array,y))

LMM(re::Vector,X::DenseMatrix,y::DataVector) = LMM(re,map(LT,re),X,convert(Array,y))

LMM(re::Vector,X::DenseMatrix,y::Vector) = LMM(re,map(LT,re),X,y)

LMM(g::PooledDataVector,y::DataVector) = LMM([ReMat(g)],y)

LMM(re::Vector,y::DataVector) = LMM(re,ones(length(y),1),y)

LMM(re::Vector,y::Vector) = LMM(re,map(LT,re),ones(length(y),1),y)

chksz(A::ReMat,λ::ParamLowerTriangular) = size(λ,1) == 1
chksz(A::VectorReMat,λ::ParamLowerTriangular) = size(λ,1) == size(A.z,1)

lowerbd(A::LMM) = mapreduce(lowerbd,vcat,A.Λ)

Base.getindex(m::LMM,s::Symbol) = mapreduce(x->x[s],vcat,m.Λ)

function Base.setindex!(m::LMM,v::Vector,s::Symbol)
    s == :θ || throw(ArgumentError("only ':θ' is meaningful for assignment"))
    lam = m.Λ
    length(v) == sum(nθ,lam) || throw(DimensionMismatch("length(v) should be $(sum(nθ,lam))"))
    A = m.A
    R = m.R
    n = size(A,1)                       # inject upper triangle of A into L
    for j in 1:n, i in 1:j
        inject!(R[i,j],A[i,j])
    end
    offset = 0
    for i in eachindex(lam)
        li = lam[i]
        nti = nθ(li)
        li[:θ] = sub(v,offset + (1:nti))
        offset += nti
        for j in i:size(R,2)
            scale!(li,R[i,j])
        end
        for ii in 1:i
            scale!(R[ii,i],li)
        end
        inflate!(R[i,i])
    end
    cfactor!(R)
end

"""
`fit(m)` -> `m`

Optimize the objective using an NLopt optimizer.
"""
function StatsBase.fit(m::LMM, verbose::Bool=false, optimizer::Symbol=:default)
    th = m[:θ]
    k = length(th)
    if optimizer == :default
        optimizer = hasgrad(m) ? :LD_MMA : :LN_BOBYQA
    end
    opt = NLopt.Opt(optimizer, k)
    NLopt.ftol_rel!(opt, 1e-12)   # relative criterion on deviance
    NLopt.ftol_abs!(opt, 1e-8)    # absolute criterion on deviance
    NLopt.xtol_abs!(opt, 1e-10)   # criterion on parameter value changes
    NLopt.lower_bounds!(opt, lowerbd(m))
    feval = 0
    geval = 0
    function obj(x::Vector{Float64}, g::Vector{Float64})
        feval += 1
        m[:θ] = x
        val = objective(m)
        if length(g) == length(x)
            geval += 1
            grad!(g,m)
        end
        val
    end
    if verbose
        function vobj(x::Vector{Float64}, g::Vector{Float64})
            feval += 1
            m[:θ] = x
            val = objective(m)
            print("f_$feval: $(round(val,5)), [")
            showcompact(x[1])
            for i in 2:length(x) print(","); showcompact(x[i]) end
            println("]")
            if length(g) == length(x)
                geval += 1
                grad!(g,m)
            end
            val
        end
        NLopt.min_objective!(opt, vobj)
    else
        NLopt.min_objective!(opt, obj)
    end
    fmin, xmin, ret = NLopt.optimize(opt, th)
    ## very small parameter values often should be set to zero
    xmin1 = copy(xmin)
    modified = false
    for i in eachindex(xmin1)
        if 0. < abs(xmin1[i]) < 1.e-5
            modified = true
            xmin1[i] = 0.
        end
    end
    if modified
        m[:θ] = xmin1
        if (ff = objective(m)) < fmin
            fmin = ff
            copy!(xmin,xmin1)
        else
            m[:θ] = xmin
        end
    end
    m.opt = OptSummary(th,xmin,fmin,feval,geval,optimizer)
    if verbose println(ret) end
    m
end

grad!(v,lmm::LMM) = v

hasgrad(lmm::LMM) = false

"""
Regenerate the last column of `m.A` from `m.trms`

This should be called after updating parts of `m.trms[end]`, typically the response.
"""
function regenerateAend!(m::LMM)
    n = Base.LinAlg.chksquare(m.A)
    trmn = m.trms[n]
    for i in 1:n
        Ac_mul_B!(m.A[i,n],m.trms[i],trmn)
    end
    m
end

"""
Reset the value of `m.θ` to the initial values
"""
function resetθ!(m::LMM)
    m[:θ] = m.opt.initial
    m.opt.feval = m.opt.geval = -1
    m
end

"Add an identity matrix to the argument, in place"
inflate!(D::Diagonal{Float64}) = (d = D.diag; for i in eachindex(d) d[i] += 1 end; D)

function inflate!{T<:FloatingPoint}(A::StridedMatrix{T})
    n = Base.LinAlg.chksquare(A)
    for i in 1:n
        @inbounds A[i,i] += 1
    end
    A
end

"""
LD(A) -> log(det(A)) for A diagonal, HBlkDiag, or UpperTriangular
"""
function LD{T}(d::Diagonal{T})
    r = log(one(T))
    dd = d.diag
    for i in eachindex(dd)
        r += log(dd[i])
    end
    r
end

function LD{T}(d::HBlkDiag{T})
    r = log(one(T))
    aa = d.arr
    p,q,k = size(aa)
    for j in 1:k, i in 1:p
        r += log(aa[i,i,j])
    end
    r
end

function LD{T}(d::DenseMatrix{T})
    r = log(one(T))
    n = Base.LinAlg.chksquare(d)
    for j in 1:n
        r += log(d[j,j])
    end
    r
end

function LD(m::LMM)
    R = m.R
    s = 0.
    for i in eachindex(m.Λ)
        s += LD(R[i,i])
    end
    2.*s
end

"Negative twice the log-likelihood"
objective(m::LMM) = LD(m) + nobs(m)*(1.+log(2π*varest(m)))

function Base.LinAlg.Ac_ldiv_B!{T}(D::UpperTriangular{T,Diagonal{T}},B::DenseMatrix{T})
    m,n = size(B)
    dd = D.data.diag
    length(dd) == m || throw(DimensionMismatch(""))
    for j in 1:n, i in 1:m
        B[i,j] /= dd[i]
    end
    B
end

function Base.LinAlg.Ac_ldiv_B!{T}(A::UpperTriangular{T,HBlkDiag{T}},B::DenseMatrix{T})
    m,n = size(B)
    aa = A.data.arr
    r,s,k = size(aa)
    m == Base.LinAlg.chksquare(A) || throw(DimensionMismatch())
    scr = Array(T,(r,n))
    for i in 1:k
        bb = sub(B,(i-1)*r+(1:r),:)
        copy!(bb,Base.LinAlg.Ac_ldiv_B!(UpperTriangular(sub(aa,:,:,i)),copy!(scr,bb)))
    end
    B
end

function Base.LinAlg.Ac_ldiv_B!{T}(D::UpperTriangular{T,Diagonal{T}},B::SparseMatrixCSC{T})
    m,n = size(B)
    dd = D.data.diag
    length(dd) == m || throw(DimensionMismatch(""))
    nzv = nonzeros(B)
    rv = rowvals(B)
    for j in 1:n, k in nzrange(B,j)
        nzv[k] /= dd[rv[k]]
    end
    B
end

Base.LinAlg.A_ldiv_B!{T<:FloatingPoint}(D::Diagonal{T},B::DenseMatrix{T}) =
    Base.LinAlg.Ac_ldiv_B!(D,B)

function Base.logdet(t::UpperTriangular)
    n = Base.LinAlg.chksquare(t)
    mapreduce(log,(+),diag(t))
end

function Base.logdet(t::LowerTriangular)
    n = Base.LinAlg.chksquare(t)
    mapreduce(log,(+),diag(t))
end

function Base.logdet{T<:FloatingPoint}(R::HBlkDiag{T})
    ra = R.arr
    ret = zero(T)
    r,s,k = size(ra)
    for i in 1:k
        ret += logdet(UpperTriangular(sub(ra,:,:,i)))
    end
    ret
end

function Base.LinAlg.A_rdiv_Bc!{T<:FloatingPoint}(A::SparseMatrixCSC{T},B::Diagonal{T})
    m,n = size(A)
    dd = B.diag
    n == length(dd) || throw(DimensionMismatch(""))
    nz = nonzeros(A)
    for j in eachindex(dd)
        @inbounds scale!(sub(nz,nzrange(A,j)),inv(dd[j]))
    end
    A
end

function Base.LinAlg.A_rdiv_Bc!{T<:FloatingPoint}(A::Matrix{T},B::Diagonal{T})
    m,n = size(A)
    dd = B.diag
    n == length(dd) || throw(DimensionMismatch(""))
    for j in eachindex(dd)
        @inbounds scale!(sub(A,:,j),inv(dd[j]))
    end
    A
end

function Base.LinAlg.A_rdiv_B!(A::StridedVecOrMat,D::Diagonal)
    m, n = size(A, 1), size(A, 2)
    if n != length(D.diag)
        throw(DimensionMismatch("diagonal matrix is $(length(D.diag)) by $(length(D.diag)) but left hand side has $n columns"))
    end
    (m == 0 || n == 0) && return A
    dd = D.diag
    for j = 1:n
        dj = dd[j]
        if dj == 0
            throw(SingularException(j))
        end
        for i = 1:m
            A[i,j] /= dj
        end
    end
    A
end

AIC(m::LMM) = objective(m) + 2npar(m)

BIC(m::LMM) = objective(m) + npar(m)*log(nobs(m))

Base.LinAlg.A_rdiv_Bc!(A::StridedVecOrMat,D::Diagonal) = A_rdiv_B!(A,D)

Base.cholfact(m::LMM) = UpperTriangular(m.R[end,end][1:end-1,1:end-1])

fixef(m::LMM) = cholfact(m)\m.R[end,end][1:end-1,end]

StatsBase.coef(m::LMM) = fixef(m)

StatsBase.nobs(m::LMM) = size(m.trms[end],1)

function StatsBase.vcov(m::LMM)
    Rinv = Base.LinAlg.inv!(cholfact(m))
    varest(m)*Rinv*Rinv'
end

"""
Number of parameters in the model.

Note that `size(m.trms[end],2)` is `length(coef(m)) + 1`, thereby accounting
for the scale parameter, σ, that is profiled out.
"""
npar(m::LMM) = size(m.trms[end],2) + length(m[:θ])

"""
returns the square root of the penalized residual sum-of-squares

This is the bottom right element of the bottom right block of m.R
"""
sqrtpwrss(m::LMM) = m.R[end,end][end,end]

"""
returns s², the estimate of σ², the variance of the conditional distribution of Y given B
"""
varest(m::LMM) = pwrss(m)/nobs(m)

"""
returns the penalized residual sum-of-squares
"""
pwrss(m::LMM) = abs2(sqrtpwrss(m))

"""
Simulate `N` response vectors from `m`, refitting the model.  The function saveresults
is called after each refit.

To save space the last column of `m.trms[end]`, which is the response vector, is overwritten
by each simulation.  The original response is restored before returning.
"""
function bootstrap(m::LMM,N::Integer,saveresults::Function,σ::Real=-1,mods::Vector{LMM}=LMM[])
    if σ < 0.
        σ = √varest(m)
    end
    trms = m.trms
    nt = length(trms)
    Xy = trms[nt]
    n,pp1 = size(Xy)
    X = sub(Xy,:,1:pp1-1)
    y = sub(Xy,:,pp1)
    y0 = copy(y)
    fev = X*coef(m)
    vfac = [σ*convert(LowerTriangular,λ) for λ in m.Λ]  # lower Cholesky factors of relative covariances
    remats = Matrix{Float64}[]
    for i in eachindex(vfac)
        vsz = size(trms[i],2)
        nr = size(vfac[i],1)
        push!(remats,reshape(zeros(vsz),(nr,div(vsz,nr))))
    end
    for kk in 1:N
        for i in 1:n
            y[i] = fev[i] + σ*randn()
        end
        for j in eachindex(remats)
            mat = vec(A_mul_B!(vfac[j],randn!(remats[j])))
            A_mul_B!(1.,trms[j],vec(mat),1.,y)
        end
        regenerateAend!(m)
        resetθ!(m)
        fit(m)
        saveresults(kk,m)
    end
    copy!(y,y0)
    regenerateAend!(m)
    resetθ!(m)
    fit(m)
end