text:
using Flux3D, Zygote, BenchmarkTools, Printf

function setup_benchmark_record(names)
    benchmarks = Dict{String,Any}()
    for name in names
        benchmarks[name] = []
    end
    return benchmarks
end

function generate_pcloud(npoints::Int)
    points = ones(Float32, 3, npoints)
    points = cumsum(points, dims = 2) / npoints
    return points
end

function generate_trimesh(npoints::Int)
    v = ones(3, npoints)
    v = cumsum(v, dims = 2) / npoints
    f = reshape(repeat(collect(1:npoints), 3), 3, npoints)
    return TriMesh([Float32.(v)], [Int32.(f)])
end

function run_benchmarks!(benchmarks, x, benchmark_func, device)
    for (func, args, name) in x
        args = args .|> device
        trial_f = @benchmark $benchmark_func($func, $args)
        back_args = [x -> func(x...), args]
        trial_t = @benchmark $benchmark_func($gradient, $back_args)
        time_f = minimum(trial_f.times) * 1.0e-6
        time_t = minimum(trial_t.times) * 1.0e-6
        time_b = time_t - time_f
        println("$name (forward): $time_f ms")
        println("$name (back): $time_b ms")
        println("$name (total): $time_t ms")
        push!(benchmarks[name], [time_f, time_b, time_t])
    end
end

npoint_arr = 2 .^ [6, 8, 10, 12, 14]

names = ["sample_points", "chamfer_distance", "edge_loss", "laplacian_loss"]

cpu_bm = setup_benchmark_record(names)
gpu_bm = setup_benchmark_record(names)

println("DEVICE: CPU")
for _npoints in npoint_arr
    arr = [
        # (sample_points, [generate_trimesh(_npoints), _npoints],
        #  "sample_points"),
        (
            chamfer_distance,
            [generate_pcloud(_npoints), generate_pcloud(_npoints)],
            "chamfer_distance",
        ),
        (edge_loss, [generate_trimesh(_npoints)], "edge_loss"),
        (laplacian_loss, [generate_trimesh(_npoints)], "laplacian_loss"),
    ]
    println("Running benchmarks for npoints = $_npoints")
    run_benchmarks!(cpu_bm, arr, (op, pc) -> op(pc...), cpu)
    println()
end

using CUDA
if has_cuda()
    println("CUDA is on. Running GPU Benchmarks")
    println("DEVICE: GPU")
    for _npoints in npoint_arr
        arr = [
            # (sample_points, [generate_trimesh(_npoints), _npoints],
            #  "sample_points"),
            (
                chamfer_distance,
                [generate_pcloud(_npoints), generate_pcloud(_npoints)],
                "chamfer_distance",
            ),
            (edge_loss, [generate_trimesh(_npoints)], "edge_loss"),
            (laplacian_loss, [generate_trimesh(_npoints)], "laplacian_loss"),
        ]
        println("Running benchmarks for npoints = $_npoints")
        run_benchmarks!(gpu_bm, arr, (op, pc) -> (CUDA.@sync op(pc...)), gpu)
        println()
    end
end

function save_bm(fname, rep, cpu_benchmarks, gpu_benchmarks)
    open(fname, "w") do io
        device = "cpu"
        for (key, values) in cpu_benchmarks
            for (p, v) in zip(npoint_arr, values)
                Printf.@printf(
                    io,
                    "Flux3D(Forward) %s %s %s %d %f ms\n",
                    rep,
                    device,
                    key,
                    p,
                    v[1]
                )
                Printf.@printf(
                    io,
                    "Flux3D(Backward) %s %s %s %d %f ms\n",
                    rep,
                    device,
                    key,
                    p,
                    v[2]
                )
                Printf.@printf(
                    io,
                    "Flux3D(Total) %s %s %s %d %f ms\n",
                    rep,
                    device,
                    key,
                    p,
                    v[3]
                )
            end
        end

        device = "gpu"
        for (key, values) in gpu_benchmarks
            for (p, v) in zip(npoint_arr, values)
                Printf.@printf(
                    io,
                    "Flux3D(Forward) %s %s %s %d %f ms\n",
                    rep,
                    device,
                    key,
                    p,
                    v[1]
                )
                Printf.@printf(
                    io,
                    "Flux3D(Backward) %s %s %s %d %f ms\n",
                    rep,
                    device,
                    key,
                    p,
                    v[2]
                )
                Printf.@printf(
                    io,
                    "Flux3D(Total) %s %s %s %d %f ms\n",
                    rep,
                    device,
                    key,
                    p,
                    v[3]
                )
            end
        end
    end
end

fname = joinpath(@__DIR__, "bm_flux3d_metrics.txt")
save_bm(fname, "Metrics", cpu_bm, gpu_bm)
@info "Benchmarks have been saved at $fname"


text:
# Term of a differential equation  coef * d^order / dt^order

struct DETerm
    coef::Union{Real, Expr, Symbol}
    order::Union{Real, Expr, Symbol}
end

# Struct to hold a list of terms for the left and right handside of the equation

struct DiffEqu
    leftvar::Symbol     # ϵ
    rightvar::Symbol    # σ
    leftde::Set{DETerm}
    rightde::Set{DETerm}
end

#  Basic constructor of an empty equation

function DiffEqu(leftvar::Symbol, rightvar::Symbol)
    return(DiffEqu(leftvar, rightvar, Set{DETerm}(), Set{DETerm}()))
end

# To add terms to the diff equation

function push!(de::DiffEqu, var::Symbol, t::DETerm)
    if var==de.leftvar
        Base.push!(de.leftde, t)
    elseif var==de.rightvar
        Base.push!(de.rightde, t)
    end
    return(nothing)
end

function push!(de::DiffEqu, var::Symbol, tt::Tuple)
    for t in tt
        push!(de,var,t)
    end
    return(nothing)
end

#  Basic constructor of an empty equation

function DiffEqu(;kwargs...)
    vars = keys(kwargs)
    de = DiffEqu(vars...)#  , values(kwargs)...)
    push!(de, vars[1], kwargs[vars[1]])
    push!(de, vars[2], kwargs[vars[2]])
    return(de)
end

# Returns a suitable derivative expression

function diffexpr(var::Symbol, order)
    if order==0
        return(var)
    elseif order==1  #(integer)
        dvar=Symbol("d"*string(var)*"dt")
        return(dvar)
    elseif typeof(order)<:Real && floor(order)==order  #(integer value)
        dvar=Symbol("d$order"*string(var)*"dt$order")
        return(:($dvar,))
    else
        return(:(dᵅdtᵅ($var, $order)))   #  May need time data too.
    end
end

# Returns a laplace/frequency power expression

function laplaceexpr(order, var=:s)
    if order==0
        return(1)
    elseif order==1  #(integer)
        return(var)
    else
        return(:($var^$order))
    end
end

#
#   Placeholder function - does nothing at the moment but will eventually return a derivative. May need time data too.
#
function dᵅ_dtᵅ(vat, order)
    return(order)
end

# Tools to operate on expressions

function addexpr(ex, tex)
    return( :($ex + $tex) )
end

function multexpr(ex, tex)
    if ex==1
        return(tex)
    elseif tex==1
        return(ex)
    else
        return( :($ex * $tex) )
    end
end

#
# Builing the differential equation
#

function builddeexpr(var::Symbol, de::Set{DETerm})
    ex = :()
    if length(de)==1
        t = first(de)
        ex = multexpr(t.coef, diffexpr(var, t.order))
    else
        s = copy(de)
        t = pop!(s)
        ex = multexpr(t.coef, diffexpr(var, t.order))
        for t in s
            ex = addexpr(ex, multexpr(t.coef, diffexpr(var, t.order)) )
        end
    end
    return(ex)
end

function builddeexpr(de::DiffEqu)
    leftex = builddeexpr(de.leftvar, de.leftde)
    rightex = builddeexpr(de.rightvar, de.rightde)
    return( :($leftex - $rightex) )
end

#
# Building the laplace / harmonic responses
#

function buildlaplaceexpr(de::Set{DETerm}, var=:s)
    ex = :()
    if length(de)==1
        t = first(de)
        ex = multexpr(t.coef, laplaceexpr(t.order, var))
    else
        s = copy(de)
        t = pop!(s)
        ex = multexpr(t.coef, laplaceexpr(t.order, var))
        for t in s
            ex = addexpr(ex, multexpr(t.coef, laplaceexpr(t.order, var)) )
        end
    end
    return(ex)
end

function builddynamicexpr(de::DiffEqu)
    leftex = buildlaplaceexpr(de.leftde, :(im * ω))
    rightex = buildlaplaceexpr(de.rightde, :(im * ω))
    return( :($leftex / $rightex) )
end

function buildrelaxationexpr(de::DiffEqu)
    leftex = buildlaplaceexpr(de.leftde)
    rightex = buildlaplaceexpr(de.rightde)
    return( :($leftex / (s * $rightex) ) )
end

function buildcreepexpr(de::DiffEqu)
    leftex = buildlaplaceexpr(de.leftde)
    rightex = buildlaplaceexpr(de.rightde)
    return( :($rightex / (s * $leftex) ) )
end

println("=====================")
println("Maxwell model")
println("=====================")

demax = DiffEqu(ϵ = DETerm(:η,1), σ = (DETerm(:k,0),DETerm(:η,1)) )

println("Differential equation (cost function)")
println(builddeexpr(demax))
println("Relaxation modulus")
println(buildrelaxationexpr(demax))
println("Creep compliance")
println(buildcreepexpr(demax))
println("Dynamic modulus")
println(builddynamicexpr(demax))

println("=====================")
println("Fract Zener model")
println("=====================")

defzener = DiffEqu( ϵ = ( DETerm(:cₐ,:α), DETerm(:cᵧ,:γ), DETerm(:(cₐ*cᵧ/cᵦ),:(α+γ-β)) ),
                    σ = (DETerm(1,0),DETerm(:(cₐ/cᵦ),:(α-β)))  )

println("Differential equation (cost function)")
println(builddeexpr(defzener))
println("Relaxation modulus")
println(buildrelaxationexpr(defzener))
println("Creep compliance")
println(buildcreepexpr(defzener))
println("Dynamic modulus")
println(builddynamicexpr(defzener))

#=

=====================
Maxwell model
=====================
Differential equation (cost function)
η * dϵdt - (η * dσdt + k * σ)
Relaxation modulus
(η * s) / (s * (η * s + k))
Creep compliance
(η * s + k) / (s * (η * s))
Dynamic modulus
(η * (im * ω)) / (η * (im * ω) + k)
=====================
Fract Zener model
=====================
Differential equation (cost function)
((((cₐ * cᵧ) / cᵦ) * dᵅdtᵅ(ϵ, (α + γ) - β) + cᵧ * dᵅdtᵅ(ϵ, γ)) + cₐ * dᵅdtᵅ(ϵ, α)) - (σ + (cₐ / cᵦ) * dᵅdtᵅ(σ, α - β))
Relaxation modulus
((((cₐ * cᵧ) / cᵦ) * s ^ ((α + γ) - β) + cᵧ * s ^ γ) + cₐ * s ^ α) / (s * (1 + (cₐ / cᵦ) * s ^ (α - β)))
Creep compliance
(1 + (cₐ / cᵦ) * s ^ (α - β)) / (s * ((((cₐ * cᵧ) / cᵦ) * s ^ ((α + γ) - β) + cᵧ * s ^ γ) + cₐ * s ^ α))
Dynamic modulus
((((cₐ * cᵧ) / cᵦ) * (im * ω) ^ ((α + γ) - β) + cᵧ * (im * ω) ^ γ) + cₐ * (im * ω) ^ α) / (1 + (cₐ / cᵦ) * (im * ω) ^ (α - β))

=#

ext: jl

max_stars_repo_name: JuliaRheology/RheoBenchAndTest

max_issues_repo_name: JuliaRheology/RheoBenchAndTest

max_forks_repo_name: JuliaRheology/RheoBenchAndTest


text:
function linearsearch(arr,x)
    n = length(arr);t=0;y=0
    for i=1:n
        if (arr[i]==x)
            t=1;y=i
        break
        end
    end
    if (t==1)
        println("element found!! and it's in position = ",y)
    else
        println("not found")
    end
end

v = [144 89 4 9 95 87 86 25]
x = 89
println("list = ",v ," and the number to be searched = ",x)

linearsearch(v,x)

max_stars_repo_name: landrumb/julia

max_issues_repo_name: landrumb/julia

max_forks_repo_name: landrumb/julia


text:
# -----------------------------------------------------------
#
# Implementing a nodal basis for DG sparse grids,
# Multidimensional case
#
# -----------------------------------------------------------

#=

    We continue the work started in 1D_Nodal_Basis.jl by applying the
    tensor product construction together with the sparse cutoff
    to constract a way to translate between the modal basis used in
    linear wave evolution and the nodal/point basis necessary for
    efficiently handling nonlinearities in the time evolution equations

=#

# Check whether (level1, cell1) overlaps with (level2, cell2)
function relevant_cell_1D(k::Int, level1::Int, cell1::Int, level2::Int, cell2::Int)
    left1  = (cell1-1)/(1<<max(0, level1-1))
    right1 = cell1/(1<<max(0, level1-1))
    left2  = (cell2-1)/(1<<max(0, level2-1))
    right2 = cell2/(1<<max(0, level2-1))
    return (left1 <= left2 && right1 >= right2) || (left1 >= left2 && right1 <= right2)
end

# Multidimensional version of the above
function relevant_cell(k::Int, level1::CartesianIndex{D},
                               cell1 ::CartesianIndex{D},
                               level2::CartesianIndex{D},
                               cell2 ::CartesianIndex{D}) where D
    bool = true
    for d in 1:D
        bool &= relevant_cell_1D(k, level1[d]-1, cell1[d], level2[d]-1, cell2[d])
    end
    return bool
end

function inner_loop(i::Int, j::Int, I::Array{Int, 1}, J::Array{Int, 1}, V::Array{Float64, 1},
                    k::Int, l2::CartesianIndex{D},
                    l1::CartesianIndex{D}, c1::CartesianIndex{D}, m1::CartesianIndex{D},
                    js_1D::CartesianIndex{D}, mat_1D::Array{T, 2}, atol::T) where {D, T<:Real}
    cells::NTuple{D, Int} = ntuple(i -> 1<<max(0, l2[i]-2), D)
    modes::NTuple{D, Int} = ntuple(i -> k, D)
    cellrange = CartesianIndices(cells)
    moderange = CartesianIndices(modes)
    for c2 in cellrange
        !relevant_cell(k, l2, c2, l1, c1) && (i += prod(modes); continue)

        for m2 in moderange
            val = one(T)
            for d in 1:D
                i_1D = get_index_1D(k, l2[d], c2[d], m2[d])
                val *= @inbounds mat_1D[i_1D, js_1D[d]]
                val == 0 && break
            end
            # (abs(val) < eps(T)) && (i += 1; continue)
            (abs(val) < atol) && (i += 1; continue)
            push!(I, i)
            push!(J, j)
            push!(V, val)
            i += 1
        end
    end
    return i
end

function make_column(j::Int, I::Array{Int, 1}, J::Array{Int, 1}, V::Array{T, 1},
                     k::Int, n::Int, l1::CartesianIndex{D}, c1::CartesianIndex{D}, m1::CartesianIndex{D},
                     mat_1D::Array{T, 2}, scheme::Val{Scheme}, atol::T) where {D, T<:Real, Scheme}
    js_1D = CartesianIndex(ntuple(d -> get_index_1D(k, l1[d], c1[d], m1[d]), D))
    levels::NTuple{D, Int} = ntuple(i -> (n+1),D)
    i = 1
    for l2 in CartesianIndices(levels)
        cutoff(scheme, l2, n) && continue
        i = inner_loop(i, j, I, J, V, k, l2, l1, c1, m1, js_1D, mat_1D, atol)
    end
end

function transform(D::Int, k::Int, n::Int, mat_1D::Array{T,2}; scheme="sparse", atol=1e-12) where {T<:Real}
    transform(Val(D), k, n, mat_1D, Val(Symbol(scheme)), atol)
end
function transform(::Val{D}, k::Int, n::Int, mat_1D::Array{T,2}, scheme::Val{Scheme}, atol) where {D, T<:Real, Scheme}
    levels::NTuple{D, Int} = ntuple(i -> (n+1),D)
    modes ::NTuple{D, Int} = ntuple(i -> k, D)
    I = Int[]
    J = Int[]
    V = T[]

    j = 1
    for l1 in CartesianIndices(levels)
        cutoff(scheme, l1, n) && continue

        cells1 = ntuple(i -> 1<<max(0, l1[i]-2), D)
        for c1 in CartesianIndices(cells1)
            for m1 in CartesianIndices(modes)
                make_column(j, I, J, V, k, n, l1, c1, m1, mat_1D, scheme, atol)
                j += 1
            end
            # Run the garbage collector relatively frequently to prevent major memory usage
            # sum(c1.I) - D % 3 == 0 && GC.gc()
        end
    end
    # return threshold(sparse(I, J, V), atol)
    return sparse(I, J, V)
end

function transform(D::Int, k::Int, n::Int, from::String, to::String; scheme="sparse", atol=1e-12)
    # It's better not to multiply the matrices from modal -> nodal with nodal -> points
    # or vice-versa. Instead, have both transform matrices in memory
    # and multiply a vector once by each of the two to get to the other basis -
    # So this if-statement is not an efficient thing to do in general
    if sort([from, to]) == ["modal", "points"]
        T1 = transform(D, k, n, from, "nodal"; scheme=scheme, atol=atol)
        T2 = transform(D, k, n, "nodal", to; scheme=scheme, atol=atol)
        return threshold(T2 * T1, atol)
    else
        mat_1D = Matrix(transform_1D(k, n, from, to))
        return transform(D, k, n, mat_1D; scheme=scheme, atol=atol)
    end
end

function make_modal2point_matrices(D::Int, k::Int, n::Int)
    println("making modal -> nodal transform (may take a while)")
    m2n = transform(D, k, n, "modal", "nodal")
    println("making nodal -> points transform")
    n2p = transform(D, k, n, "nodal", "points")
    return m2n, n2p
end

function make_point2modal_matrices(D::Int, k::Int, n::Int)
    println("making points -> nodal transform")
    p2n = transform(D, k, n, "points", "nodal")
    println("making nodal -> modal transform (may take a while)")
    n2m = transform(D, k, n, "nodal", "modal")
    return p2n, n2m
end

# --------------------------------------------------
# Below this line is for comparison testing only -
# Not for use
# --------------------------------------------------

function transform2(D::Int, k::Int, n::Int, mat_1D::Array{T,2}; scheme="sparse", atol=1e-12) where {T<:Real}
    transform2(Val(D), k, n, mat_1D, Val(scheme), atol)
end
function transform2(::Val{D}, k::Int, n::Int, mat_1D::Array{T,2}, scheme::Val{Scheme}, atol) where {D, T<:Real, Scheme}
    levels = ntuple(i->(n+1),D)
    modes  = ntuple(i -> k, D)
    hier_ref = D2Vref(1, k, n)
    I = Int[]; J = Int[]; V = Float64[]
    i::Int = 1; j::Int = 1; val::Float64 = one(Float64)
    index::Int = 1

    levelrange = CartesianIndices(levels); moderange = CartesianIndices(modes)
    js_1D::Array{Int,1} = zeros(Int, D)
    cells1::Array{Int,1} = zeros(Int, D);
    cells2::Array{Int,1} = zeros(Int, D)
    for l1 in levelrange
        cutoff(scheme, l1, n) && continue

        cells1 = [1<<max(0, l1[q]-2) for q in 1:D]
        for c1 in CartesianIndices((cells1...,))
            for m1 in moderange
                #1D indices corresponding to each (l1[d], c1[d], m1[d])
                js_1D = [get_index_1D(k, l1[d], c1[d], m1[d]) for d in 1:D]

                i = 1
                for l2 in levelrange
                    cutoff(schemel, l2, n) && continue

                    cells2 = [1<<max(0, l2[q]-2) for q in 1:D]
                    for c2 in CartesianIndices((cells2...,))
                        # If there's no overlap - skip
                        !relevant_cell(k, l2, c2, l1, c1) && (i += prod(modes); continue)
                        for m2 in moderange
                            val = one(Float64)
                            for d in 1:D
                                i_1D = get_index_1D(k, l2[d], c2[d], m2[d])
                                val *= mat_1D[i_1D, js_1D[d]]
                                val == 0 && break
                            end
                            # (abs(val) < eps(T)) && (i += 1; continue)
                            (abs(val) < atol) && (i += 1; continue)
                            push!(I, i)
                            push!(J, j)
                            push!(V, val)
                            i += 1
                        end
                    end
                end
                j += 1
            end
        end
    end
    # return threshold(sparse(I, J, V), atol)
    return sparse(I, J, V)
end

function transform2(D::Int, k::Int, n::Int, from::String, to::String; scheme="sparse", atol=1e-12)
    if Set{String}([from, to]) == Set{String}(["modal", "points"])
        T1 = transform2(D, k, n, from, "nodal"; scheme=scheme, atol=atol)
        T2 = transform2(D, k, n, "nodal", to; scheme=scheme, atol=atol)
        return threshold(T2 * T1, atol)
    else
        mat_1D = Matrix(transform_1D(k, n, from, to))
        return transform2(D, k, n, mat_1D; scheme=scheme, atol=atol)
    end
end


text:
# Practical Machine learning
# Support Vector Machine example
# Chapter 6

## Commands from https://gist.github.com/2833165
## gcc -Wall -fPIC -c *.c
## cc -shared -Wl,-soname,libtest.so -o libtest.so *.o

## Adapt to svm:
## gcc -Wall -fPIC -c *.c *.cpp
## cc -shared -Wl,-soname,libsvm.so -o libsvm.so *.o

libsvm = dlopen("reference/libsvm.so")

function svm_nodes(X::Array{Float64,2})  ## TODO:  write version for sparse X
    r = size(X)[1]
    c = size(X)[2]
    print(r)
    print(c)
    x = reshape(transpose(X), (r * c, 1))
    svn = ccall(dlsym(libsvm, :sparsify),
                     Ptr{Void},
                     (Ptr{Float64}, Int32, Int32),
                     x, r, c)
    return svn
end

function svm_problem(y::Array{Float64,1}, X::Array{Float64,2})
    svn = svm_nodes(X)
    l = int32(size(y)[1])
    svmprob = ccall(dlsym(libsvm, :_jl_svm_problem),
                    Ptr{Void},
                    (Ptr{Float64}, Int32, Ptr{Void}),
                    y, l, svn)
    return svmprob
end
function svm_problem(f::Formula, df::DataFrame)
    mm = model_matrix(f, df)
    y = mm.response[:,1]
    X = mm.model
    svmprob = svm_problem(y, X)
    return svmprob
end
svm_problem(f::Expr, df::DataFrame) = svm_problem(Formula(f), df)

const svm_type_table = {"C_SVC"=>1,
           "NU_SVC"=>2,
           "ONE_CLASS"=>3,
           "EPSILON_SVR"=>4,
           "NU_SVR"=>5}
const kern_type_table = {"LINEAR"=>1,
            "POLY"=>2,
            "RBF"=>3,
            "SIGMOID"=>4,
            "PRECOMPUTED"=>5}

function svmparameter(svm_type::String, kernel_type::String, degree::Int32,
                      gamma::Float64, coef0::Float64, cache::Float64, tolerance::Float64,
                      cost::Float64, nu::Float64, #nweights,
                      p::Float64, shrinking::Int32, probability::Int32)
    svm_type_int = get(svm_type_table, uppercase(svm_type), 0)
    kernel_type_int = get(kern_type_table, uppercase(svm_type), 0)
    svmpar = ccall(dlsym(libsvm, :_jl_svm_par),
                   Ptr{Void},
                   (Int32, Int32, Int32,
                    Float64, Float64, Float64, Float64,
                    Float64, Float64,
                    Float64, Int32, Int32),
                   svm_type_int, kernel_type_int, degree,
                   gamma, coef0, cache, tolerance,
                   cost, nu, #nweights,
                   p, shrinking, probability)
    return svmpar
end

function svmtrain(prob, param)
    svm_model = ccall(dlsym(libsvm, :svm_train),
                      Ptr{Void},
                      (Ptr{Void}, Ptr{Void}),
                      prob, param)
end

function svmpredict(model, X)
    svn = svmnode(X)
    r = size(X)[1]
    res = Array(Float64, (r,))
    ccall(dlsym(libsvm, :_jl_svm_predict),
          Ptr{Float64},
          (Ptr{Void}, Ptr{Void}, Int32, Ptr{Float64}),
          model, svn, r, res)
    return res
end


text:
function compute_int_J(params, α)
    (; r, D, W, dim) = params
    m = length(D)
    R = r[m]

    λ = α^2

    if α > 1e-6
        V, _ = compute_v(params, α, 0, R)
        J = W[m] * V / λ / R
    else
        J = 1 / dim
    end

    J
end


text:
putAdj(A::Assoc,AA::AbstractSparseMatrix) =  Assoc(copy(A.row),copy(A.col),copy(A.val),AA)
putAdj(A::Assoc,AA::Array{Int,2}) =  putAdj(A::Assoc,sparse(AA))

putRow(A::Assoc,ARow::Array{Union{AbstractString,Number},1}) =  Assoc(ARow,copy(A.col),copy(A.val),copy(A.A))

putCol(A::Assoc,ACol::Array{Union{AbstractString,Number},1}) =  Assoc(copy(A.row),ACol,copy(A.val),copy(A.A))

putVal(A::Assoc,AVal::Array{Union{AbstractString,Number},1}) =  Assoc(copy(A.row),copy(A.col),AVal,copy(A.A))


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

# Test integer conversion routines from int.jl

for y in (-4, Float32(-4), -4.0, big(-4.0))
    @test flipsign(3, y)  == -3
    @test flipsign(-3, y) == 3
    @test copysign(3, y)  == -3
    @test copysign(-3, y) == -3
end

for y in (4, Float32(4), 4.0, big(4.0))
    @test flipsign(3, y)  == 3
    @test flipsign(-3, y) == -3
    @test copysign(3, y)  == 3
    @test copysign(-3, y) == 3
end

# Result type must be type of first argument
for T in (Base.BitInteger_types..., BigInt,
          Rational{Int}, Rational{BigInt},
          Float16, Float32, Float64)
    for U in (Base.BitInteger_types..., BigInt,
              Rational{Int}, Rational{BigInt},
              Float16, Float32, Float64)
        @test typeof(copysign(T(3), U(4))) === T
        @test typeof(flipsign(T(3), U(4))) === T
    end
end

for s1 in (-1,+1), s2 in (-1,+1)
    @test flipsign(Int16(3s1), Float16(3s2)) === Int16(3s1*s2)
    @test flipsign(Int32(3s1), Float32(3s2)) === Int32(3s1*s2)
    @test flipsign(Int64(3s1), Float64(3s2)) === Int64(3s1*s2)
end

@test signed(3) == 3
@test signed(UInt(3)) == 3
@test isa(signed(UInt(3)), Int)
@test signed(UInt(0) - 1) == -1
@test_throws InexactError signed(UInt(-3))

@test bswap(Int8(3)) == 3
@test bswap(UInt8(3)) == 3
@test bswap(Int16(3)) == 256*3
@test bswap(Int16(256)) == 1
@test bswap(Int16(257)) == 257
@test bswap(Int32(1)) == 2^(3*8)
@test bswap(Int32(2)^(3*8)) == 1
@test bswap(Int64(1)) == Int64(2)^(7*8)
@test bswap(Int64(2)^(7*8)) == 1
@test bswap(Int128(1)) == Int128(2)^(15*8)
@test bswap(Int128(2)^(15*8)) == Int128(1)
@test bswap(UInt128(2)^(15*8)) == UInt128(1)

@test count_zeros(10) == Sys.WORD_SIZE - 2
@test count_zeros(UInt8(10)) == 6

@test convert(Signed, UInt128(3)) === Int128(3)
@test convert(Signed, false) === 0
@test convert(Signed, true) === 1

for (II, UU) in ((Int8,UInt8), (Int16,UInt16),
                 (Int32,UInt32), (Int64,UInt64),
                 (Int128,UInt128))
    @test convert(Unsigned, II(3)) === UU(3)
end

for T in (Float32, Float64)
    @test convert(Unsigned, T(3.0)) === UInt(3)
end

@test trunc(3) == 3
@test trunc(Integer, 3) == 3

@test floor(3) == 3
@test ceil(3) == 3

@test big"2"^100 == BigInt(2)^100
@test isa(big"2", BigInt)
@test big"1.0" == BigFloat(1.0)
@test_throws ArgumentError big"1.0.3"
@test_throws ArgumentError big"pi"

@test round(UInt8, 123) == 123
@test mod(123, UInt8) == 0x7b

primitive type MyBitsType <: Integer 8 end
@test_throws MethodError ~reinterpret(MyBitsType, 0x7b)

UItypes = Base.BitUnsigned_types
SItypes = Base.BitSigned_types

for T in UItypes, S in UItypes
    @test promote(S(3), T(3)) === (sizeof(T) < sizeof(S) ? (S(3), S(3)) : (T(3), T(3)))
end

for T in SItypes, S in SItypes
    @test promote(S(3), T(3)) === (sizeof(T) < sizeof(S) ? (S(3), S(3)) : (T(3), T(3)))
end

for T in SItypes, S in UItypes
    R = sizeof(S) < sizeof(Int) ? Int : S
    @test promote(R(3), T(3)) === (sizeof(R) < sizeof(T) ? (T(3), T(3)) : (R(3), R(3)))
end

# Test limiting conversions
for T in (Int8, Int16, Int32, Int64)
    max_val = Int128(typemax(T))
    @test convert(T, max_val) == max_val
    @test_throws InexactError convert(T, max_val+1)

    m = Int128(typemin(T))
    @test convert(T, m) == m
    @test_throws InexactError convert(T, m-1)
end

for T in (UInt8, UInt16, UInt32, UInt64)
    max_val = Int128(typemax(T))
    @test convert(T, max_val) == max_val
    @test_throws InexactError convert(T, max_val+1)
    @test_throws InexactError convert(T, -1)
end

# Test bit shifts
for T in Base.BitInteger_types
    nbits = 8*sizeof(T)
    issigned = typemin(T) < 0
    highbit = T(2) ^ (nbits-1)
    val1 = 0x1234567890abcdef % T % highbit
    val2 = val1 + highbit
    for val in (val1, val2)
        for count in 0:nbits+1
            ucount, scount = unsigned(count), signed(count)
            # Note: We assume modulo semantics for the arithmetic operations
            # used here
            if count < nbits
                @test val << ucount === val * T(2)^count
                @test val >>> ucount ===
                    fld(unsigned(val), unsigned(T(2))^count) % T
            else
                @test val << ucount === T(0)
                @test val >>> ucount === T(0)
            end
            @test val << scount === val << ucount
            @test val << -scount === val >> ucount
            @test val >>> scount === val >>> ucount
            @test val >>> -scount === val << ucount
            if count < (issigned ? nbits-1 : nbits)
                @test val >> ucount === fld(val, T(2)^count)
            else
                @test val >> ucount === T(val<0 ? -1 : 0)
            end
            @test val >> scount === val >> ucount
            @test val >> -scount === val << ucount
        end
    end
end

@test widen(UInt8(3)) === UInt32(3)
@test widen(UInt16(3)) === UInt32(3)
@test widen(UInt32(3)) === UInt64(3)
@test widen(UInt64(3)) === UInt128(3)
@test widen(UInt128(3)) == 3
@test typeof(widen(UInt128(3))) == BigInt

@test widen(Int8(-3)) === Int32(-3)
@test widen(Int16(-3)) === Int32(-3)
@test widen(Int32(-3)) === Int64(-3)
@test widen(Int64(-3)) === Int128(-3)
@test widen(Int128(-3)) == -3
@test typeof(widen(Int128(-3))) == BigInt

@test widemul(false, false) == false
@test widemul(false, 3) == 0
@test widemul(3, true) == widemul(true, 3) == 3

let i=Int64(2)^63-1, k=widemul(i,i)
    @test widemul(i,i)==85070591730234615847396907784232501249
    j=div(k,2)
    @test div(k,j)==2
    j=div(k,5)
    @test rem(k,j)==4
end

# issue #3596
@test Int128(1)<<0 == 1
@test repr(Int128(1)<<1) == "2"

# issue #16700
@test_throws MethodError 1.0 >> 8

# PR #16988
@test true << 2 === 1 << 2
@test true >> 2 === 1 >> 2
@test true >>> 2 === 1 >>> 2

@test @inferred(unsafe_trunc(Int8, 127)) === Int8(127)
@test unsafe_trunc(Int8, 128) === Int8(-128)
@test unsafe_trunc(Int8, -127) === Int8(-127)
@test unsafe_trunc(Int8, -128) === Int8(-128)
@test unsafe_trunc(Int8, -129) === Int8(127)


text:
#=
Create all the colormaps from cmocean and register them with
Matplotlib.

Example:
	include("path/to/cmocean.jl")

	imshow(Z, cmap="ice")

WARNING: This script overwrites the "gray" colormap from Matplotlib.
         The one from cmocean has a slightly longer range. If you
         don't want that, just delete "gray" below.

=#
using PyPlot
using CSV

for name in ["algae", "amp", "balance", "curl", "deep", "delta",
	"dense", "gray", "haline", "ice", "matter", "oxy",
	"phase", "solar", "speed", "tempo", "thermal", "turbid"]

	path = "$(@__DIR__)/rgb/$name-rgb.txt"
	df = CSV.read(path,header=false,delim=' ')
	raw = convert(Array{Float64,2},df)

	cm = ColorMap(name, raw)
	plt[:register_cmap](name=name,cmap=cm)

	cm = ColorMap(name*"_r", raw[end:-1:1,:])
	plt[:register_cmap](name=name*"_r",cmap=cm)
end


text:
import Mercury as Hg
import Mercury: LoopRateLimiter
using LinearAlgebra
using StaticArrays
using Statistics
using Test

# @testset "Rate Limited Loops" begin
function runrate(f, lrl::LoopRateLimiter, time_s = 1)
    N = round(Int, lrl.rate * time_s)
    times = zeros(N)
    Hg.reset!(lrl)
    for i = 1:N
        times[i] = @elapsed begin
            Hg.startloop(lrl)
            f()
            sleep(lrl)
        end
    end
    return 1 ./ times
end

@generated function runrate_macro1(f, ::Val{rate}, time_s = 1) where {rate}
    quote
        N = round(Int, $rate * time_s)
        lrl = LoopRateLimiter(rate)
        t_start = time_ns()
        Hg.@rate for i = 1:N
            f()
        end $rate
        t_elapsed = (time_ns() - t_start) * 1e-9
        return N / t_elapsed
    end
end

function runrate_macro2(f, rate, time_s = 1)
    N = round(Int, rate * time_s)
    lrl = LoopRateLimiter(rate)
    t_start = time_ns()
    Hg.@rate for i = 1:N
        f()
    end lrl
    t_elapsed = (time_ns() - t_start) * 1e-9
    return N / t_elapsed
end

@generated function runrate_macro3(f, ::Val{rate}, time_s = 1) where {rate}
    quote
        N = round(Int, $rate * time_s)
        lrl = LoopRateLimiter(rate)
        t_start = time_ns()
        i = 0
        Hg.@rate while (i < N)
            f()
            i += 1
        end $rate
        t_elapsed = (time_ns() - t_start) * 1e-9
        return N / t_elapsed
    end
end

function runrate_macro4(f, rate, time_s = 1)
    N = round(Int, rate * time_s)
    lrl = LoopRateLimiter(rate)
    t_start = time_ns()
    i = 0
    Hg.@rate while (i < N)
        f()
        i += 1
    end lrl
    t_elapsed = (time_ns() - t_start) * 1e-9
    return N / t_elapsed
end

##
function mykernel()
    A = @SMatrix randn(10, 10)
    B = @SMatrix randn(10, 10)
    C = A'A + B'B
    x = @SVector randn(10)
    cholA = cholesky(Symmetric(C))
    C \ x
end
1 / @elapsed mykernel()

@testset "Rate limiter tests" begin
    rate = 1000
    test_time = 0.1 # sec
    lrl = LoopRateLimiter(rate)
    @test Hg.doesblock(lrl)
    error = (rate - median(runrate(mykernel, lrl, test_time))) / rate * 100
    @test error < 10  # less than 5% error

    error = (rate - runrate_macro1(mykernel, Val(rate), test_time)) / rate * 100
    @test error < 10  # less than 5% error
    error = (rate - runrate_macro2(mykernel, rate, test_time)) / rate * 100
    @test error < 10  # less than 5% error
    error = (rate - runrate_macro3(mykernel, Val(rate), test_time)) / rate * 100
    @test error < 10  # less than 5% error
    error = (rate - runrate_macro4(mykernel, rate, test_time)) / rate * 100
    @test error < 10  # less than 5% error

    rate = 40
    test_time = 0.1 # sec
    thresh = 6  # percent error
    lrl = LoopRateLimiter(rate)
    @test !Hg.doesblock(lrl)
    runrate(mykernel, lrl, test_time)
    error = (rate - median(runrate(mykernel, lrl, test_time))) / rate * 100
    @test error < thresh

    error = (rate - runrate_macro1(mykernel, Val(rate), test_time)) / rate * 100
    @test error < thresh
    error = (rate - runrate_macro2(mykernel, rate, test_time)) / rate * 100
    @test error < thresh
    error = (rate - runrate_macro3(mykernel, Val(rate), test_time)) / rate * 100
    @test error < thresh
    error = (rate - runrate_macro4(mykernel, rate, test_time)) / rate * 100
    @test error < thresh

end


text:
##codegen.jl
# functions relating to verilog code generation and verilation.

function normalized_key(mod::Union{Function, Symbol}, p::Tuple)
  m_rep = Symbol(mod)

  for k in keys(__global_dependency_cache)
    if (k[1] == Symbol(mod)) && ([p...] == [pair.second for pair in k[2:end]])
      return k
    end
  end

  throw(ArgumentError("Key not found!"))
end

doc"""
  `Verilog.find_dependencies(mod, (params...))`

  outputs the set of dependencies of a certain module.
"""
@generated function find_dependencies(mod, p::Tuple, exclude = [])
  if mod <: Function  #then it's a "manual call"
    tuplegen = :(key = normalized_key(mod, p))
  elseif mod == Symbol
    tuplegen = :(key = (mod, p...))
  else
    throw(ArgumentError("module type unsupported: $mod"))
  end

  quote
    $tuplegen
    deps = __global_dependency_cache[key]

    result = Tuple[]
    #now iterate through the dependencies.

    for d in deps
      if !((d in result) || (d in exclude))
        unshift!(result, d)
        result = union(result, find_dependencies(d[1], d[2:end]))
      end
    end
    result
  end
end

const lintoff_stmts = """
/* verilator lint_off DECLFILENAME */
/* verilator lint_off UNUSED */
/* verilator lint_off PINMISSING */

"""

doc"""
  `generate_verilog_file(path, mod, (params...); exclude = [])`

  pass this function a path, (verilog module) name and the parameters, and it
  puts that module and its dependencies into a file in the desired path.  The
  filename will be (module name + suffix).v, with module name & suffix as
  generated by the interal specifications.
"""
function generate_verilog_file(path::String, mod::Function, p::Tuple; exclude = [])
  #test to see if the path exists.
  isdir(path) || mkdir(path)
  #first call the function and make it.
  modulestring = mod(p...)
  modulename = ""

  #next call its dependencies.
  deps = find_dependencies(mod, p, exclude)
  depstrings = [__global_definition_cache[dep].txt for dep in deps]

  #next, retrieve the proper name for the file.
  nkey = normalized_key(mod, p)

  modulename = __global_definition_cache[nkey].module_name

  #next, create a file for the header.
  vfile_fio = open(string(path, "/", modulename, ".v"), "w")
  try
    #write the contents
    println(vfile_fio, lintoff_stmts)
    println(vfile_fio, join(depstrings, "\n\n"))
    println(vfile_fio, modulestring)
  finally
    #close the directory
    close(vfile_fio)
  end

  #return its dependencies so that we can keep track of all modules we've
  #created so far.
  unshift!(deps, nkey)
  deps
end

export generate_verilog_file

doc"""
  `verilate(mod::Function, p::Tuple)`
  uses verilator to create a .so file corresponding to the library.
  Each of the functions in the mods array will be the top-level element
  in a freestanding .v file.  They may refer to each other.
  the "mods" are a tuple of the function and the parameter.

  Passed options:
    path::String      - defaults to "." - the library .so file will be placed in the path.
    libname::String   - file name for the library.  Defaults to libVerilated-[function_with_suffix].so
    with_source::Bool - defaults to "false" - should the c source be included?
                        this is useful if you're building a library around the
                        verilated module.
"""
function verilate(mod::Function, p::Tuple; path::String = ".", libname::String="", with_source = false)
  #check to see if the path_to_c_library is actually a dir.
  isdir(path) || mkdir(path)

  nkey = normalized_key(mod, p)
  #retrieve the proper name of this function
  name = __global_definition_cache[nkey].module_name

  if libname == ""
    libname = "libVerilated-$(name).so"
  end

  println("verilating to generate $path/$libname")

  #create the temp directory
  mktempdir((tdir)->begin
    #create the interface cpp file.
    open((file) -> begin
      println(file, verilator_adapter(nkey))
    end, string(tdir, "/", name, "_interface.cpp"), "w")

    println("tdir:", tdir)
    generate_verilog_file(tdir, mod, p)

    #next move into the temporary directory and then
    cd(()->begin
      run(`verilator -Wall --cc -CFLAGS -fPIC $name.v --top-module $name --exe $(name)_interface.cpp`)
    end, tdir)

    cd(()-> begin
      #run the make file.  This will throw an error due to not finding a --main
      #but that's OK it just needs to make the object files.
      run(`make -f V$name.mk`)

      object_files = ["$tdir/obj_dir/$f" for f in readdir() if f[end-1:end] == ".o"]

      #manually link the object files into a shared library.
      run(`cc -Wall -shared $object_files -o $libname`)
    end, "$tdir/obj_dir")

    #in the case that the caller requests the source code
    if (with_source)
      #copy all the source files to the destination directory
      cp("$tdir/", "$path/"; remove_destination=true)
    else
      #finally, save the temporary verilated file into the main path, squashing it if necessary.
      mv("$tdir/obj_dir/$libname", "$path/$libname"; remove_destination=true)
    end
  end)
end

doc"""
  `@verilate function parameters [path]`

  verilates a @verilog function with the passed parameter and places the shared
  object library into [path].  By default, the path used is going to be
  "./libVerilated-[function].so".

  the function is then hooked up to a function with signature []
"""
macro verilate(f, p = (), path...)
  #f should be a symbol
  @assert isa(f, Symbol) "@verilate must be run on a function symbol"
  #check to make sure that f is in the function database.

  @assert isa(eval(Main, f), Function) "@verilate must be run on a function symbol"

  if length(path) == 1
    path = path[1]
    @assert isa(path, String) "@verilate paths must be strings"
  else
    path = "./libVerilated-$f.so"
  end

  (dir, fname) = splitdir(path)
  realdir = realpath(dir)

  #trigger compilation of the function verilog by executing the parameter tuple
  #on a blank version of the function.
  v_fn = eval(Main, f)
  v_p  = eval(Main, p)
  v_fn(v_p...)

  #assign the module cache
  mod_cache = Verilog.__global_definition_cache[normalized_key(f, v_p)]
  name = mod_cache.module_name

  #generate the input tuple and input signature
  input_signature = :(())
  input_converts = :()
  input_ccall = :(ccall((:set, $path), Void, $input_signature))
  for idx = 1:length(mod_cache.inputs)
    iname = Symbol(:input, idx)
    input_converts = :($input_converts; $iname = UInt64($iname))
    push!(input_signature.args, :UInt64)
    push!(input_ccall.args, iname)
  end

  #create "c" versions of the function
  cfunctionname = Symbol(name, "_c")

  #build the output type object to be consistent with the c definition.
  #also build the return tuple command.
  otypename = Symbol("__", name, "_output_type")
  output_type_generator = :(type $otypename; end)
  return_tuple = :(())
  create_tuple = :(())
  for idx = 1:length(mod_cache.outputlist)
    oname = Symbol(:output, idx)
    push!(output_type_generator.args[3].args, :($oname::UInt64))
    push!(return_tuple.args, :(output_ref.$oname))
    push!(create_tuple.args, 0x0000_0000_0000_0000)
  end

  #generate the code for the evaluation function
  if length(mod_cache.outputlist) == 1
    efunc = :(function ($cfunctionname)()
      #create an output_ref
      $input_ccall
      ccall((:step, $path), Void, ())
      res = ccall((:get, $path), UInt64, ())
      res
    end)
  else
    efunc = :(function ($cfunctionname)()
      #create an output_ref
      output_ref = $otypename(($create_tuple)...)
      $input_ccall
      ccall((:step, $path), Void, ())
      ccall((:get, $path), Void, (Ptr{$otypename},), pointer_from_objref(output_ref))
      $return_tuple
    end)
  end

  #append input definitions onto the efunc.
  for idx = 1:length(mod_cache.inputs)
    input_name = Symbol(:input, idx)
    input_stmt = :($input_name::Unsigned)
    push!(efunc.args[1].args, input_stmt)
  end

  #this needs to be eval'd because julia is finicky.
  eval(Main, output_type_generator)

  #this also needs to be eval'd because of the way that julia
  #parses and executes ccalls.
  eval(Main, :(verilate($f, $p, path = $realdir, libname = $fname)))

  esc(quote
    #initialize the verilated version of the function.
    ccall((:init, $path), Void, ())

    $efunc
  end)
end

export verilate
export @verilate


text:
module TestMacrosWhere

import Test: @testset, @test
import Volcanito: @where, materialize
import DataFrames: DataFrame

@testset "@where" begin
    df = DataFrame(
        a = [1, 2, 3, 4],
        b = [2, 3, 4, 5],
    )

    plan1 = @where(df, a < b)
    res1 = materialize(plan1)
    @test res1 == df

    plan2 = @where(df, a > 1)
    res2 = materialize(plan2)
    @test res2 == df[2:4, :]
end

end


text:
function process_likelihood_model(ρ_list::Vector,Eₘ_list::Vector)
    # Generate the A matrix used to calculate likelihoods
    # The A matrix depends on the input states and measurement operators
    dimsmatch(ρ_list,Eₘ_list)
    sum(abs,data(sum(Eₘ_list))-I)<1E-15 ||
        throw(ArgumentError("Eₘ operators do not form a valid POVM!"))
    prep = (Eₘ,ρ) -> transpose(vec(Array(ρ⊗transpose(Eₘ))))
    return reduce(vcat,vec(prep.(Eₘ_list,permutedims(ρ_list))))
end

# below is the actual projected gradient descent algorithm from Knee, et al.
# Quantum process tomography via completely positive and trace-preserving
# projection. Phys. Rev. A 98, 062336 (2018).

function pgd_process_tomo(M::Matrix{T}, A::Matrix; tol=1E-10, cptp_tol=1E-8, info=false, cbfun=nothing) where T
    # Choi process matrix reconstruction by maximum likelihood projected gradient descent
    size(A,1)==length(M) || throw(ArgumentError("A matrix inconsistent with number of measurements!"))
    abs(sum(M)-1)<1/4 || throw(ArgumentError("measurement counts not normalized!"))
    # infer space dimensions from A matrix
    d = isqrt(isqrt(size(A,2)))
    # initial Choi matrix guess, the identity map
    C = Matrix{Complex{T}}(I/d,d^2,d^2)
    # objective and gradient functions setup
    f = C -> loglikelihood(M,C,A)
    ∇f = C -> loglikelihood_gradient(M,C,A)
    # metaparameters & initial cost calculation
    μ⁻¹ = 2d^2/T(3); γ = T(0.3)
    c₁ = T(1E6); c₂ = f(C)
    info && println("start cost = $c₂")
    # iterate through projected gradient descent steps, with backtracking
    h = CPTP_helpers(C)
    stp = 0
    while c₁ - c₂ > tol
        stp += 1
        c₁, ∇C = c₂, ∇f(C)
        D = project_CPTP(C .- μ⁻¹.*∇C, h, cptp_tol) - C
        α = one(T); Π = γ*real(D⋅∇C)
        while (c₂ = f(C .+ α.*D)) > c₁ + α*Π
            α = α/2 # backtrack
        end
        C .= C .+ α.*D
        if cbfun !== nothing # run callback function; e.g. to calc fidelity at each step
            cbfun(stp,c₂,C)
        end
    end
    info && println("final cost = $c₂, Δ = $(c₁-c₂), number of steps: $stp")
    return C
end

function loglikelihood(M::Matrix{T}, C::Matrix{Complex{T}}, A::Matrix) where {T<:Real}
    # Binomial statistics for the measurement count probability, up to some irrelevant constant
    logP = log.(max.(real.(A*vec(C)), T(1E-12)))
    return -(M⋅logP)
end

function loglikelihood_gradient(M::Matrix{T}, C::Matrix{Complex{T}}, A::Matrix) where {T<:Real}
    mMP = .-vec(M)./max.(real.(A*vec(C)), T(1E-12))
    return unvec(A'*mMP)
end

function project_CPTP(C::Matrix{Complex{T}}, h::Tuple, tol::Real=1E-8) where {T<:Real}
    # generate storage objects
    X₁ = copy(C); Y₁ = zero(X₁);
    X₂ = copy(Y₁); Y₂ = copy(Y₁)
    P = copy(Y₁); Q = copy(Y₁)
    ΔP = one(T); ΔQ = one(T)
    D, V, MdagM, 𝕀 = h
    # iterate through TP & CP projections
    while ΔP^2 + ΔQ^2 + 2*abs(P⋅X₂-P⋅X₁) + 2*abs(Q⋅Y₂-Q⋅Y₁) > tol
        project_TP!(Y₂, X₁+P, MdagM, 𝕀)
        ΔP = norm2_diff(X₁,Y₂)
        P .= X₁ .- Y₂ .+ P
        project_CP!(X₂, Y₂+Q, D, V)
        ΔQ = norm2_diff(Y₂,X₂)
        Q .= Y₂ .- X₂ .+ Q
        X₁, X₂ = X₂, X₁
        Y₁, Y₂ = Y₂, Y₁
    end
    return X₁
end

function project_CP!(X::Matrix{Complex{T}}, C::Matrix{Complex{T}}, D::Vector{T}, V::Matrix) where {T<:Real}
    # Project the process onto the completely positive subspace by making the
    # Choi matrix positive semidefinite
    # We do this by taking the eigendecomposition, setting any negative
    # eigenvalues to 0, and reconstructing the Choi matrix
    hermfact!(D,V,Hermitian(hermitianize!(C)))
    D .= max.(D,zero(T))
    mul!(C,Diagonal(D),V')
    mul!(X,V,C)
end

function project_TP!(Y::Matrix{Complex{T}}, C::Matrix{Complex{T}}, MdagM::Matrix, 𝕀::Matrix) where {T<:Real}
    # Project the process onto the trace-preserving subspace
    d⁻¹ = one(T)/isqrt(size(C,1))
    mul!(vec(Y), MdagM, vec(C))
    Y .= C .- d⁻¹.*(Y.-𝕀)
end

function CPTP_helpers(C::Matrix{Complex{T}}) where {T<:Real}
    D = Vector{T}(undef,size(C,1))
    V = Matrix{Complex{T}}(undef,size(C))
    MdagM, 𝕀 = TP_helper_matrices(C)
    return D, V, MdagM, 𝕀
end

function TP_helper_matrices(C::Matrix{Complex{T}}) where {T<:Real}
    d = isqrt(size(C,1))
    𝕀 = Matrix{T}(I,d,d); k = zeros(Complex{T},1,d)
    # this could be done more efficiently, but it doesn't matter
    M = sum(i->(k[i]=1; k[mod1(i-1,d)]=0; (𝕀 ⊗ k) ⊗ (𝕀 ⊗ k)), 1:d)
    return M'*M, 𝕀⊗𝕀
end


text:
cd(@__DIR__); allNGWPs = false; include("setups/toronto.jl")
gr(dpi = 200)

## fdensity 16 most important VM-NGWP vectors (ignore the DC vector)
VM_NGWP = vm_ngwp(𝚽, GP_dual)
f = zeros(N); for i in 1:N; f[i] = length(findall(dist_X[:,i] .< 1/minimum(edge_weight))); end
G_Sig.f = reshape(f, (N, 1))
dmatrix_VM = ngwp_analysis(G_Sig, VM_NGWP)
dvec_vm_ngwp, BS_vm_ngwp = ngwp_bestbasis(dmatrix_VM, GP_dual)
important_idx = sortperm(dvec_vm_ngwp[:].^2; rev = true)
println("================fdensity-VM-NGWP-top-basis-vectors=================")
for i in 2:17
    dr, dc = BS_vm_ngwp.levlist[important_idx[i]]
    w = VM_NGWP[dr, dc, :]
    j, k, l = NGWP_jkl(GP_dual, dr, dc)
    print("(j, k, l) = ($(j), $(k), $(l))  ")
    if j == jmax
        print("φ_{$(GP_dual.ind[dr]-1)}")
    end
    println()
    gplot(A, X; width=1)
    scatter_gplot!(X; marker = w, plotOrder = :s2l, ms = 3)
    plt = plot!(cbar = false, clims = (-0.075,0.075))
    savefig(plt,
        "../figs/Toronto_fdensity_DAG_VM_NGWP_ibv$(lpad(i,2,"0"))_j$(j)_k$(k)_l$(l).png")
end


text:
module OpenStreetMapXPlot

using OpenStreetMapX

import Plots
import PyPlot

export plotmap, addroute!, plot_nodes! #Plotting

include("plot.jl") #plotting
include("layers.jl") #colors

end #module


text:
using Distributions, Gadfly, Random
using Turing, MCMCChains
Turing.turnprogress(false)  # turn off the progress monitor

Random.seed!(3)

# Construct the data points
N = 30
μs = [[0., 0.], [5., 5.]]
x1 = rand(MvNormal(μs[1], 1.), N)
x2 = rand(MvNormal(μs[2], 1.), N)
x = hcat([x1, x2]...)
plot(x=x[1,:], y=x[2,:])

# Define the probabilistic model
@model GaussianMixtureModel(x) = begin
    K, N = size(x)

    # Global variables
    μ = 0.
    σ = 1.
    μ1 ~ Normal(μ, σ)
    μ2 ~ Normal(μ, σ)
    locs = [μ1, μ2]
    α = 1.
    weights ~ Dirichlet(K, α)  # latent
    # weights = ones(K) / K;  # observed

    # Local context
    z = Vector{Int64}(undef, N)
    for i in 1:N
        z[i] ~ Categorical(weights)
        x[:,i] ~ MvNormal([locs[z[i]], locs[z[i]]], 1.)
    end
    return z  # only return *latent* local variables
end

gmm_model = GaussianMixtureModel(x);

# Perform MCMC inference
gmm_sampler = Gibbs(
    PG(100, :z),  # Particle Gibbs (discrete)
    HMC(0.05, 10, :μ1, :μ2, :weights)  # Hamiltonian Monte Carlo (continuous)
    # HMC(0.05, 10, :μ1, :μ2, )  # Hamiltonian Monte Carlo (continuous)
);
@time chain = sample(gmm_model, gmm_sampler, 100);

ids = findall(map(name -> occursin("μ", name), names(chain)));
p=plot(chain[:, ids, :], legend=true, labels = ["Mu 1" "Mu 2"], colordim=:parameter)

chain = chain[:, :, 1];

function predict(x, y, w, μ)
    # Use log-sum-exp trick for numeric stability.
    return Turing.logaddexp(
        log(w[1]) + logpdf(MvNormal([μ[1], μ[1]], 1.), [x, y]),
        log(w[2]) + logpdf(MvNormal([μ[2], μ[2]], 1.), [x, y])
    )
end

contour(range(-5, stop = 3), range(-6, stop = 2),
    (x, y) -> predict(x, y, [0.5, 0.5], [mean(chain[:μ1].value), mean(chain[:μ2].value)])
)
scatter!(x[1,:], x[2,:], legend = false, title = "Synthetic Dataset")

assignments = collect(skipmissing(mean(chain[:k].value, dims=1).data))
scatter(x[1,:], x[2,:],
    legend = false,
    title = "Assignments on Synthetic Dataset",
    zcolor = assignments)


text:
export Dirichlet, Neumann

function boundary_state!(
    nf,
    land::LandModel,
    state⁺::Vars,
    diff⁺::Vars,
    aux⁺::Vars,
    n̂,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    bctype,
    t,
    _...,
)
    args = (state⁺, diff⁺, aux⁺, n̂, state⁻, diff⁻, aux⁻, bctype, t)
    soil_boundary_state!(nf, land, land.soil, land.soil.water, args...)
    soil_boundary_state!(nf, land, land.soil, land.soil.heat, args...)
end

function boundary_state!(
    nf,
    land::LandModel,
    state⁺::Vars,
    aux⁺::Vars,
    n̂,
    state⁻::Vars,
    aux⁻::Vars,
    bctype,
    t,
    _...,
)
    args = (state⁺, aux⁺, n̂, state⁻, aux⁻, bctype, t)
    soil_boundary_state!(nf, land, land.soil, land.soil.water, args...)
    soil_boundary_state!(nf, land, land.soil, land.soil.heat, args...)
end

abstract type AbstractBoundaryFunctions end

"""
    struct Dirichlet{Fs, Fb} <: AbstractBoundaryFunctions

A concrete type to hold the surface state and bottom state variable
values/functions, if Dirichlet boundary conditions are desired.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct Dirichlet{Fs, Fb} <: AbstractBoundaryFunctions
    "Surface state boundary condition"
    surface_state::Fs = nothing
    "Bottom state boundary condition"
    bottom_state::Fb = nothing
end

"""
    struct Neumann{Fs, Fb} <: AbstractBoundaryFunctions

A concrete type to hold the surface and/or bottom diffusive flux
values/functions, if Neumann boundary conditions are desired.

Note that these are intended to be scalar values. In the boundary_state!
functions, they are multiplied by the `ẑ` vector (i.e. the normal vector `n̂`
to the domain at the upper boundary, and -`n̂` at the lower boundary. These
normal vectors point out of the domain.)

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct Neumann{Fs, Fb} <: AbstractBoundaryFunctions
    "Surface flux boundary condition"
    surface_flux::Fs = nothing
    "Bottom flux boundary condition"
    bottom_flux::Fb = nothing
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
# Functional Interface
_issues_dir() = joinpath(bookdir(), "Issues")

_issue_file(refid::String) = joinpath(_read_reports_dir(), string(refid, ".txt"))


text:
mutable struct TigerPOMDP <: POMDP{Bool, Int64, Bool}
    r_listen::Float64
    r_findtiger::Float64
    r_escapetiger::Float64
    p_listen_correctly::Float64
    discount_factor::Float64
end
TigerPOMDP() = TigerPOMDP(-1.0, -100.0, 10.0, 0.85, 0.95)

states(::TigerPOMDP) = (true, false)
observations(::TigerPOMDP) = (true, false)

stateindex(::TigerPOMDP, s::Bool) = Int64(s) + 1
actionindex(::TigerPOMDP, a::Int) = a + 1
obsindex(::TigerPOMDP, o::Bool) = Int64(o) + 1

initial_belief(::TigerPOMDP) = DiscreteBelief(2)

const TIGER_LISTEN = 0
const TIGER_OPEN_LEFT = 1
const TIGER_OPEN_RIGHT = 2

const TIGER_LEFT = true
const TIGER_RIGHT = false

# Resets the problem after opening door; does nothing after listening
function transition(pomdp::TigerPOMDP, s::Bool, a::Int64)
    p = 1.0
    if a == 1 || a == 2
        p = 0.5
    elseif s
        p = 1.0
    else
        p = 0.0
    end
    return BoolDistribution(p)
end

function observation(pomdp::TigerPOMDP, a::Int64, sp::Bool)
    pc = pomdp.p_listen_correctly
    p = 1.0
    if a == 0
        sp ? (p = pc) : (p = 1.0-pc)
    else
        p = 0.5
    end
    return BoolDistribution(p)
end

function observation(pomdp::TigerPOMDP, s::Bool, a::Int64, sp::Bool)
    return observation(pomdp, a, sp)
end

function reward(pomdp::TigerPOMDP, s::Bool, a::Int64)
    r = 0.0
    a == 0 ? (r+=pomdp.r_listen) : (nothing)
    if a == 1
        s ? (r += pomdp.r_findtiger) : (r += pomdp.r_escapetiger)
    end
    if a == 2
        s ? (r += pomdp.r_escapetiger) : (r += pomdp.r_findtiger)
    end
    return r
end
reward(pomdp::TigerPOMDP, s::Bool, a::Int64, sp::Bool) = reward(pomdp, s, a)

initialstate_distribution(pomdp::TigerPOMDP) = BoolDistribution(0.5)

actions(::TigerPOMDP) = [0,1,2]

function upperbound(pomdp::TigerPOMDP, s::Bool)
    return pomdp.r_escapetiger
end

discount(pomdp::TigerPOMDP) = pomdp.discount_factor

function initialobs(p::TigerPOMDP, s::Bool, rng::AbstractRNG)
    d = observation(p, 0, s) # listen
    return rand(rng, d)
end


text:
export alloc_special, alloc_local

"Allocates on-device memory statically from the specified address space."
@generated function alloc_special(::Val{id}, ::Type{T}, ::Val{as}, ::Val{len}) where {id,T,as,len}
    Context() do ctx
        eltyp = convert(LLVMType, T; ctx)

        # old versions of GPUArrays invoke _shmem with an integer id; make sure those are unique
        if !isa(id, String) || !isa(id, Symbol)
            id = "alloc_special_$id"
        end

        T_ptr_i8 = convert(LLVMType, LLVMPtr{T,as}; ctx)

        # create a function
        llvm_f, _ = create_function(T_ptr_i8)

        # create the global variable
        mod = LLVM.parent(llvm_f)
        gv_typ = LLVM.ArrayType(eltyp, len)
        gv = GlobalVariable(mod, gv_typ, string(id), as)
        if len > 0
            linkage!(gv, LLVM.API.LLVMExternalLinkage)
            # NOTE: Backend doesn't support initializer for local AS
        end

        # by requesting a larger-than-datatype alignment, we might be able to vectorize.
        # TODO: Make the alignment configurable
        alignment!(gv, Base.max(32, Base.datatype_alignment(T)))

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            ptr_with_as = gep!(builder, gv, [ConstantInt(0; ctx),
                                             ConstantInt(0; ctx)])
            ptr = bitcast!(builder, ptr_with_as, T_ptr_i8)
            ret!(builder, ptr)
        end

        call_function(llvm_f, LLVMPtr{T,as})
    end
end

@inline alloc_local(id, T, len) = alloc_special(Val(id), T, Val(AS.Local), Val(len))

@inline @generated function alloc_string(::Val{str}) where str
    Context() do ctx
        T_pint8 = LLVM.PointerType(LLVM.Int8Type(ctx), AS.Global)
        llvm_f, _ = create_function(T_pint8)
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)
            str_ptr = globalstring_ptr!(builder, String(str))
            ptr = addrspacecast!(builder, str_ptr, T_pint8)
            ret!(builder, ptr)
        end
        call_function(llvm_f, LLVMPtr{UInt8,AS.Global})
    end
end

# TODO: Support various types of len
@inline @generated function memcpy!(dest_ptr::LLVMPtr{UInt8,DestAS}, src_ptr::LLVMPtr{UInt8,SrcAS}, len::LT) where {DestAS,SrcAS,LT<:Union{Int64,UInt64}}
    Context() do ctx
        T_nothing = LLVM.VoidType(ctx)
        T_pint8_dest = convert(LLVMType, dest_ptr; ctx)
        T_pint8_src = convert(LLVMType, src_ptr; ctx)
        T_int64 = convert(LLVMType, len; ctx)
        T_int1 = LLVM.Int1Type(ctx)

        llvm_f, _ = create_function(T_nothing, [T_pint8_dest, T_pint8_src, T_int64])
        mod = LLVM.parent(llvm_f)
        T_intr = LLVM.FunctionType(T_nothing, [T_pint8_dest, T_pint8_src, T_int64, T_int1])
        intr = LLVM.Function(mod, "llvm.memcpy.p$(DestAS)i8.p$(SrcAS)i8.i64", T_intr)
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            dest_ptr_i8 = parameters(llvm_f)[1]
            src_ptr_i8  = parameters(llvm_f)[2]
            call!(builder, intr, [dest_ptr_i8, src_ptr_i8, parameters(llvm_f)[3], ConstantInt(T_int1, 0)])
            ret!(builder)
        end
        call_function(llvm_f, Nothing, Tuple{LLVMPtr{UInt8,DestAS},LLVMPtr{UInt8,SrcAS},LT}, :dest_ptr, :src_ptr, :len)
    end
end
memcpy!(dest_ptr::LLVMPtr{T,DestAS}, src_ptr::LLVMPtr{T,SrcAS}, len::Integer) where {T,DestAS,SrcAS} =
    memcpy!(reinterpret(LLVMPtr{UInt8,DestAS}, dest_ptr), reinterpret(LLVMPtr{UInt8,SrcAS}, src_ptr), UInt64(len))
@inline @generated function memset!(dest_ptr::LLVMPtr{UInt8,DestAS}, value::UInt8, len::LT) where {DestAS,LT<:Union{Int64,UInt64}}
    Context() do ctx
        T_nothing = LLVM.VoidType(ctx)
        T_pint8_dest = convert(LLVMType, dest_ptr; ctx)
        T_int8 = convert(LLVMType, value; ctx)
        T_int64 = convert(LLVMType, len; ctx)
        T_int1 = LLVM.Int1Type(ctx)

        llvm_f, _ = create_function(T_nothing, [T_pint8_dest, T_int8, T_int64])
        mod = LLVM.parent(llvm_f)
        T_intr = LLVM.FunctionType(T_nothing, [T_pint8_dest, T_int8, T_int64, T_int1])
        intr = LLVM.Function(mod, "llvm.memset.p$(DestAS)i8.i64", T_intr)
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            call!(builder, intr, [parameters(llvm_f)[1], parameters(llvm_f)[2], parameters(llvm_f)[3], ConstantInt(T_int1, 0)])
            ret!(builder)
        end
        call_function(llvm_f, Nothing, Tuple{LLVMPtr{UInt8,DestAS},UInt8,LT}, :dest_ptr, :value, :len)
    end
end
memset!(dest_ptr::LLVMPtr{T,DestAS}, value::UInt8, len::Integer) where {T,DestAS} =
    memset!(convert(LLVMPtr{UInt8,DestAS}, dest_ptr), value, UInt64(len))


text:
using YaoBase, YaoArrayRegister
import StaticArrays: SMatrix
export RotationGate, Rx, Ry, Rz, rot

"""
    RotationGate{N, T, GT <: AbstractBlock{N, Complex{T}}} <: PrimitiveBlock{N, Complex{T}}

RotationGate, with GT both hermitian and isreflexive.

# Definition

```math
\\mathbf{I} cos(θ / 2) - im sin(θ / 2) * mat(U)
```
"""
mutable struct RotationGate{N, T <: Real, GT <: AbstractBlock{N}} <: PrimitiveBlock{N}
    block::GT
    theta::T
    function RotationGate{N, T, GT}(block::GT, theta) where {N, T, GT <: AbstractBlock{N}}
        #ishermitian(block) && isreflexive(block) ||:1
        #    throw(ArgumentError("Gate type $GT is not hermitian or not isreflexive." * " " * string(ishermitian(block)) * " " * string(isreflexive(block)) ))
        new{N, T, GT}(block, T(theta))
    end
end

RotationGate(block::GT, theta::T) where {N, T <: Real, GT<:AbstractBlock{N}} = RotationGate{N, T, GT}(block, theta)

# convert to float if theta is not a floating point
#RotationGate(block::AbstractBlock, theta) = RotationGate(block, Float64(theta))

# bindings
"""
    Rx(theta)

Return a [`RotationGate`](@ref) on X axis.

# Example

```jldoctest
julia> Rx(0.1)
rot(X gate, 0.1)
```
"""
Rx(theta) = RotationGate(X, theta)

"""
    Ry(theta)

Return a [`RotationGate`](@ref) on Y axis.

# Example

```jldoctest
julia> Ry(0.1)
rot(Y gate, 0.1)
```
"""
Ry(theta) = RotationGate(Y, theta)

"""
    Rz(theta)

Return a [`RotationGate`](@ref) on Z axis.

# Example

```jldoctest
julia> Rz(0.1)
rot(Z gate, 0.1)
```
"""
Rz(theta) = RotationGate(Z, theta)

"""
    rot(U, theta)

Return a [`RotationGate`](@ref) on U axis.
"""
rot(axis::AbstractBlock, theta) = RotationGate(axis, theta)

content(x::RotationGate) = x.block
# General definition
function mat(::Type{T}, R::RotationGate{N}) where {N, T}
    I = IMatrix{1<<N, T}()
    return I * cos(R.theta / 2) - im * sin(R.theta / 2) * mat(T, R.block)
    #return I * cos(T(R.theta) / 2) - im * sin(T(R.theta) / 2) * mat(T, R.block)
end

# Specialized
mat(::Type{T}, R::RotationGate{1, <:Any, <:XGate}) where T =
    T[cos(R.theta/2) -im * sin(R.theta/2); -im * sin(R.theta/2) cos(R.theta/2)]
mat(::Type{T}, R::RotationGate{1, <:Any, <:YGate}) where T =
     T[cos(R.theta/2) -sin(R.theta/2); sin(R.theta/2) cos(R.theta/2)]
# mat(R::RotationGate{1, T, ZGate{Complex{T}}}) where T =
#     SMatrix{2, 2, Complex{T}}(cos(R.theta/2)-im*sin(R.theta/2), 0, 0, cos(R.theta/2)+im*sin(R.theta/2))

function apply!(r::ArrayReg, rb::RotationGate)
    v0 = copy(r.state)
    apply!(r, rb.block)
    # NOTE: we should not change register's memory address,
    # or batch operations may fail
    r.state .= -im*Complex{typeof(rb.theta)}(sin(rb.theta/2))*r.state + Complex{typeof(rb.theta)}(cos(rb.theta/2))*v0
    #r.state .= -im*(sin(rb.theta/2))*r.state + (cos(rb.theta/2))*v0
    return r
end

# parametric interface
niparams(::Type{<:RotationGate}) = 1
getiparams(x::RotationGate) = x.theta
setiparams!(r::RotationGate, param::Real) = (r.theta = param; r)

YaoBase.isunitary(r::RotationGate) = true

Base.adjoint(blk::RotationGate) = RotationGate(blk.block, -blk.theta)
Base.copy(R::RotationGate) = RotationGate(R.block, R.theta)
Base.:(==)(lhs::RotationGate{TA, GTA}, rhs::RotationGate{TB, GTB}) where {TA, TB, GTA, GTB} = false
Base.:(==)(lhs::RotationGate{TA, GT}, rhs::RotationGate{TB, GT}) where {TA, TB, GT} = lhs.theta == rhs.theta

cache_key(R::RotationGate) = R.theta


text:
import Base.length
import Base.ndims

type Document
    wordVector::Array{Int64,1}
end

type TopicArray
    Assigned::Array{Int64,1}
end

type SparseCounts
    nonzeroDocIndex::Array{Int64,1}
end

type SparseDocument
    wordVector::Array{Int64,1}
    tokenIndex::Array{Int64,1}
    tokenCount::Array{Int64,1}
    docLength::Int64
end

type sparseSCVBO_params_type
    Alpha::Array{Float64,2}
    Beta::Array{Float64,2}
    BetaSum::Float64
    numTopics::Int64
    numWords::Int64
    wallClockTime::Float64
    topicCounts::Array{Float64,2}
    wordTopicCounts::Array{Float64,2}
    numCachedSamples::Int64
    cachedSamples::Array{Int64,2}
    cached_qw::Array{Float64,2}
    sampleIndx::Array{Int64,2}
    cached_Qw::Array{Float64,2}
end

sparseSCVBO_params_type(Alpha, Beta, numTopics,numWords,numDocumentIterations,numDocuments) = sparseSCVBO_params_type(Alpha,Beta,sum(Beta),
                            numTopics,numWords,0.0,ones(1,numTopics),ones(numWords,numTopics),numTopics,zeros(Int64,numWords,numTopics),
                            zeros(numWords,numTopics),ones(Int64,numWords,1),zeros(numWords,1))
length(x::sparseSCVBO_params_type) = 1
ref(x::sparseSCVBO_params_type,i::Int) = x
ndims(x::sparseSCVBO_params_type) = 0


text:
###########################################################
#
#   Generic matrix generation functions
#
###########################################################

function matrix_from_adjpairs!(a::AbstractMatrix, g::AbstractGraph, gen)
    @graph_requires g vertex_list vertex_map

    if implements_edge_list(g)
        if is_directed(g)
            for e in edges(g)
                u = source(e, g)
                v = target(e, g)
                ui = vertex_index(u, g)
                vi = vertex_index(v, g)
                a[ui, vi] = get(gen, g, u, v)
            end
        else
            for e in edges(g)
                u = source(e, g)
                v = target(e, g)
                ui = vertex_index(u, g)
                vi = vertex_index(v, g)
                val = get(gen, g, u, v)
                a[ui, vi] = val
                a[vi, ui] = val
            end
        end

    elseif implements_adjacency_list(g)
        for u in vertices(g)
            ui = vertex_index(u, g)
            for v in out_neighbors(u, g)
                vi = vertex_index(v, g)
                val = get(gen, g, u, v)
                a[ui, vi] = val
            end
        end
    else
        error("g does not implement required interface.")
    end

    return a
end

matrix_from_adjpairs(g::AbstractGraph, gen) =
    (n = num_vertices(g); matrix_from_adjpairs!(zeros(eltype(gen), n, n), g, gen))

function sparse_matrix_from_adjpairs(g::AbstractGraph, gen)
    @graph_requires g vertex_list vertex_map

    n = num_vertices(g)
    m = num_edges(g)
    ne = is_directed(g) ? m : 2m
    I = Array(Int, ne)
    J = Array(Int, ne)
    vals = Array(eltype(gen), ne)
    idx = 0

    if implements_edge_list(g)
        if is_directed(g)
            for e in edges(g)
                u = source(e, g)
                v = target(e, g)
                ui = vertex_index(u, g)
                vi = vertex_index(v, g)
                val = get(gen, g, u, v)

                idx += 1
                I[idx] = ui
                J[idx] = vi
                vals[idx] = val
            end
        else
            for e in edges(g)
                u = source(e, g)
                v = target(e, g)
                ui = vertex_index(u, g)
                vi = vertex_index(v, g)
                val = get(gen, g, u, v)

                idx += 1
                I[idx] = ui
                J[idx] = vi
                vals[idx] = val

                idx += 1
                I[idx] = vi
                J[idx] = ui
                vals[idx] = val
            end
        end

    elseif implements_adjacency_list(g)
        for u in vertices(g)
            ui = vertex_index(u, g)
            for v in out_neighbors(u, g)
                vi = vertex_index(v, g)
                val = get(gen, g, u, v)

                idx += 1
                I[idx] = ui
                J[idx] = vi
                vals[idx] = val
            end
        end
    end

    @assert(idx == ne)
    sparse(I, J, vals, n, n)
end

function matrix_from_edges!(a::AbstractMatrix, g::AbstractGraph, gen)
    @graph_requires g vertex_list vertex_map

    if implements_edge_list(g)
        if is_directed(g)
            for e in edges(g)
                u = source(e, g)
                v = target(e, g)
                ui = vertex_index(u, g)
                vi = vertex_index(v, g)
                val = get(gen, g, e)
                a[ui, vi] = val
            end
        else
            for e in edges(g)
                u = source(e, g)
                v = target(e, g)
                ui = vertex_index(u, g)
                vi = vertex_index(v, g)
                val = get(gen, g, e)
                a[ui, vi] = val
                a[vi, ui] = val
            end
        end

    elseif implements_incidence_list(g)
        for u in vertices(g)
            ui = vertex_index(u, g)
            for e in out_edges(u, g)
                v = target(e, g)
                vi = vertex_index(v, g)
                val = get(gen, g, e)
                a[ui, vi] = val
            end
        end
    else
        error("g does not implement required interface.")
    end

    return a
end

matrix_from_edges(g::AbstractGraph, gen) =
    (n = num_vertices(g); matrix_from_edges!(zeros(eltype(gen), n, n), g, gen))

function sparse_matrix_from_edges(g::AbstractGraph, gen)
    @graph_requires g vertex_list vertex_map

    n = num_vertices(g)
    m = num_edges(g)
    ne = is_directed(g) ? m : 2m
    I = Array(Int, ne)
    J = Array(Int, ne)
    vals = Array(eltype(gen), ne)
    idx = 0

    if implements_edge_list(g)
        if is_directed(g)
            for e in edges(g)
                u = source(e, g)
                v = target(e, g)
                ui = vertex_index(u, g)
                vi = vertex_index(v, g)
                val = get(gen, g, e)

                idx += 1
                I[idx] = ui
                J[idx] = vi
                vals[idx] = val
            end
        else
            for e in edges(g)
                u = source(e, g)
                v = target(e, g)
                ui = vertex_index(u, g)
                vi = vertex_index(v, g)
                val = get(gen, g, e)

                idx += 1
                I[idx] = ui
                J[idx] = vi
                vals[idx] = val

                idx += 1
                I[idx] = vi
                J[idx] = ui
                vals[idx] = val
            end
        end

    elseif implements_incidence_list(g)
        for u in vertices(g)
            ui = vertex_index(u, g)
            for e in out_edges(u, g)
                v = target(e, g)
                vi = vertex_index(v, g)
                val = get(gen, g, e)

                idx += 1
                I[idx] = ui
                J[idx] = vi
                vals[idx] = val
            end
        end
    end

    @assert(idx == ne)
    sparse(I, J, vals, n, n)
end

###########################################################
#
#   Specific matrices from graphs
#
###########################################################

### adjacency matrix

type _GenUnit{T} end

Base.get{T,V}(::_GenUnit{T}, g::AbstractGraph{V}, u::V, v::V) = one(T)
Base.eltype{T}(::_GenUnit{T}) = T

adjacency_matrix{T<:Number}(g::AbstractGraph, ::Type{T}) = matrix_from_adjpairs(g, _GenUnit{T}())
adjacency_matrix(g::AbstractGraph) = adjacency_matrix(g, Bool)

adjacency_matrix_sparse{T<:Number}(g::AbstractGraph, ::Type{T}) = sparse_matrix_from_adjpairs(g, _GenUnit{T}())
adjacency_matrix_sparse(g::AbstractGraph) = adjacency_matrix_sparse(g, Bool)

### weight matrix

type _GenEdgeWeight{Weights}
    weights::Weights
end

_GenEdgeWeight(weights::AbstractVector) = _GenEdgeWeight{typeof(weights)}(weights)

Base.get{V,E}(gen::_GenEdgeWeight, g::AbstractGraph{V,E}, e::E) = gen.weights[edge_index(e, g)]
Base.eltype(gen::_GenEdgeWeight) = eltype(gen.weights)

weight_matrix(g::AbstractGraph, weights::AbstractVector) = matrix_from_edges(g, _GenEdgeWeight(weights))
weight_matrix_sparse(g::AbstractGraph, weights::AbstractVector) = sparse_matrix_from_edges(g, _GenEdgeWeight(weights))

### distance matrix

function init_distancemat{T<:Number}(n::Int, dinf::T)
    a = fill(dinf, n, n)
    @inbounds for i = 1:n
        a[i,i] = zero(T)
    end
    return a
end

distance_matrix(g::AbstractGraph, dists::AbstractVector, dinf) =
    (n = num_vertices(g); matrix_from_edges!(init_distancemat(n, dinf), g, _GenEdgeWeight(dists)))

distance_matrix{T<:Real}(g::AbstractGraph, dists::AbstractVector{T}) =
    distance_matrix(g, dists, typemax(T))

###########################################################
#
#   Laplacian matrix
#
###########################################################

laplacian_matrix(g::AbstractGraph) = laplacian_matrix(g, Float64)

function laplacian_matrix{T<:Number}(g::AbstractGraph, ::Type{T})
    @graph_requires g vertex_list vertex_map
    !is_directed(g) || error("g must be undirected.")

    n = num_vertices(g)
    a = zeros(T, n, n)

    if implements_edge_list(g)
        for e in edges(g)
            u = source(e, g)
            v = target(e, g)
            ui = vertex_index(u, g)
            vi = vertex_index(v, g)
            if ui != vi
                a[ui, vi] -= 1
                a[vi, ui] -= 1
                a[ui, ui] += 1
                a[vi, vi] += 1
            end
        end
    elseif implements_adjacency_list(g)
        for u in vertices(g)
            ui = vertex_index(u, g)
            for v in out_neighbors(u, g)
                vi = vertex_index(v, g)
                if ui < vi
                    a[ui, vi] -= 1
                    a[vi, ui] -= 1
                    a[ui, ui] += 1
                    a[vi, vi] += 1
                end
            end
        end
    else
        error("g does not implement proper interface.")
    end
    return a
end

function laplacian_matrix{T<:Number}(g::AbstractGraph, eweights::AbstractVector{T})
    @graph_requires g vertex_list vertex_map edge_map
    !is_directed(g) || error("g must be undirected.")

    n = num_vertices(g)
    a = zeros(T, n, n)

    if implements_edge_list(g)
        for e in edges(g)
            u = source(e, g)
            v = target(e, g)
            ui = vertex_index(u, g)
            vi = vertex_index(v, g)
            if ui != vi
                wi = eweights[edge_index(e, g)]
                a[ui, vi] -= wi
                a[vi, ui] -= wi
                a[ui, ui] += wi
                a[vi, vi] += wi
            end
        end
    else
        for u in vertices(g)
            ui = vertex_index(u, g)
            for e in out_edges(u, g)
                v = target(e, g)
                vi = vertex_index(v, g)
                if ui < vi
                    wi = eweights[edge_index(e, g)]
                    a[ui, vi] -= wi
                    a[vi, ui] -= wi
                    a[ui, ui] += wi
                    a[vi, vi] += wi
                end
            end
        end
    end
    return a
end

laplacian_matrix_sparse(g::AbstractGraph) = laplacian_matrix_sparse(g, Float64)

function laplacian_matrix_sparse{T<:Number}(g::AbstractGraph, ::Type{T})
    @graph_requires g vertex_list vertex_map
    !is_directed(g) || error("g must be undirected.")

    n = num_vertices(g)
    a = zeros(T, n, n)

    nnz = num_edges(g) * 2 + n
    I = Array(Int, nnz)
    J = Array(Int, nnz)
    vals = Array(T, nnz)
    idx = 0

    degs = zeros(T, n)

    if implements_edge_list(g)
        for e in edges(g)
            u = source(e, g)
            v = target(e, g)
            ui = vertex_index(u, g)
            vi = vertex_index(v, g)
            if ui != vi
                idx += 1
                I[idx] = ui
                J[idx] = vi
                vals[idx] = -1

                idx += 1
                I[idx] = vi
                J[idx] = ui
                vals[idx] = -1

                degs[ui] += 1
                degs[vi] += 1
            end
        end
    elseif implements_adjacency_list(g)
        for u in vertices(g)
            ui = vertex_index(u, g)
            for v in out_neighbors(u, g)
                vi = vertex_index(v, g)
                if ui < vi
                    idx += 1
                    I[idx] = ui
                    J[idx] = vi
                    vals[idx] = -1

                    idx += 1
                    I[idx] = vi
                    J[idx] = ui
                    vals[idx] = -1

                    degs[ui] += 1
                    degs[vi] += 1
                end
            end
        end
    else
        error("g does not implement proper interface.")
    end

    for i = 1:n
        di = degs[i]
        if isnz(di)
            idx += 1
            I[idx] = i
            J[idx] = i
            vals[idx] = di
        end
    end

    if idx < nnz
        I = I[1:idx]
        J = J[1:idx]
        vals = vals[1:idx]
    end
    sparse(I, J, vals, n, n)
end

function laplacian_matrix_sparse{T<:Number}(g::AbstractGraph, eweights::AbstractVector{T})
    @graph_requires g vertex_list vertex_map edge_map
    !is_directed(g) || error("g must be undirected.")

    n = num_vertices(g)
    a = zeros(T, n, n)

    nnz = num_edges(g) * 2 + n
    I = Array(Int, nnz)
    J = Array(Int, nnz)
    vals = Array(T, nnz)
    idx = 0

    degs = zeros(T, n)

    if implements_edge_list(g)
        for e in edges(g)
            u = source(e, g)
            v = target(e, g)
            ui = vertex_index(u, g)
            vi = vertex_index(v, g)
            if ui != vi
                wi = eweights[edge_index(e, g)]
                idx += 1
                I[idx] = ui
                J[idx] = vi
                vals[idx] = -wi

                idx += 1
                I[idx] = vi
                J[idx] = ui
                vals[idx] = -wi

                degs[ui] += wi
                degs[vi] += wi
            end
        end
    elseif implements_adjacency_list(g)
        for u in vertices(g)
            ui = vertex_index(u, g)
            for e in out_edges(u, g)
                v = target(e, g)
                vi = vertex_index(v, g)
                if ui < vi
                    wi = eweights[edge_index(e, g)]
                    idx += 1
                    I[idx] = ui
                    J[idx] = vi
                    vals[idx] = -wi

                    idx += 1
                    I[idx] = vi
                    J[idx] = ui
                    vals[idx] = -wi

                    degs[ui] += wi
                    degs[vi] += wi
                end
            end
        end
    else
        error("g does not implement proper interface.")
    end

    for i = 1:n
        di = degs[i]
        if isnz(di)
            idx += 1
            I[idx] = i
            J[idx] = i
            vals[idx] = di
        end
    end

    if idx < nnz
        I = I[1:idx]
        J = J[1:idx]
        vals = vals[1:idx]
    end
    sparse(I, J, vals, n, n)
end

###########################################################
#
#   adjacency list from matrix
#
###########################################################

function sparse2adjacencylist{Tv,Ti<:Integer}(A::SparseMatrixCSC{Tv,Ti})
    colptr = A.colptr
    rowval = A.rowval
    n = size(A, 1)
    adjlist = Array(Array{Ti,1}, n)
    s = 0
    for j in 1:n
        adjj = Ti[]
        sizehint(adjj, colptr[j+1] - colptr[j] - 1)
        for k in colptr[j]:(colptr[j+1] - 1)
            rvk = A.rowval[k]
            if rvk != j push!(adjj, rvk) end
        end
        s += length(adjj)
        adjlist[j] = adjj
    end
    GenericAdjacencyList{Ti, Range1{Ti}, Vector{Vector{Ti}}}(!ishermitian(A),
                                                             one(Ti):convert(Ti,n),
                                                             s, adjlist)
end


text:
"""
Examples:
```julia
d"egen :x = mean(:y)"
d"egen :x = :y .+ :z"
d"egen :x = :y  if :z .> 1.0 "
d"bysort mygroup (myindex): egen :x = mean(:y)  if :z .> 1.0 "
```
`egen` is faster when combined with type declarations, e.g.
```julia
d"bysort mygroup (myindex): egen :x::Float64 = mean(:y)  if :z .> 1.0 "
```
If the type of the new column is not declared like that, `egen` will try to infer it by applying the right-hand side of the assignment
operation to the full DataFrame and taking the type of the resulting column. All columns that result from `egen` support `Missing`s.

As in Stata, `if` conditions restrict the set of rows that are being used for the RHS *and* the set of observations
that are being assigned to.

Stata's `egen` is inconsistent on how missing values are handled. Some `egen` commands (like `mean()`) ignore all rows
that have missing values for one of the variables used on the right-hand side. Other commands (like `rowtotal()`)
by default use all rows.
The default behavior of Douglass's `egen` is to ignore all rows that have `missing` entries for
one or more variables in the right-hand side expression. That does _not_ include columns used in the `if` statement.

Differences to Stata:
  - Type declarations are done via `d"egen :x::Float64 = :y .+ :z"`. Types are from Julia, of course (which means we can finally use 64-bit Integers!)
"""
macro egenerate(t::Symbol,
    by::Union{Vector{Symbol}, Nothing},
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Union{Vector{Symbol},Union{Expr, Nothing}},
    filter::Union{Expr, Nothing},
    use::Union{String, Nothing},
    options::Union{Dict{String,Any}, Nothing})
    error("""\n
        The syntax is:
            egenerate <var> = <expression> [if <expression>]
        or
            bysort <varlist> (<varlist>): egenerate <var> = <expression> [if <expression>]
    """)
end

# short form. This is a copy of the above for `gen`.
macro egen(t::Symbol,
    by::Union{Vector{Symbol}, Nothing},
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Union{Vector{Symbol},Union{Expr, Nothing}},
    filter::Union{Expr, Nothing},
    use::Union{String, Nothing},
    options::Union{Dict{String,Any}, Nothing})
    return esc(
        quote
            Douglass.@egenerate($t, $by, $sort, $arguments, $filter, $use, $options)
        end
    )
end

macro egenerate(t::Symbol,
    by::Vector{Symbol},
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Expr,
    filter::Union{Expr, Nothing},
    use::Nothing,
    options::Union{Dict{String,Any}, Nothing})
    return esc(
        quote
            Douglass.@generate_byvec!($t, $by, $sort, $arguments, $filter, $options )
        end
    )
end
# version without `by` but with `if`
macro egenerate(t::Symbol,
    by::Nothing,
    sort::Union{Vector{Symbol}, Nothing},
    arguments::Expr,
    filter::Union{Expr, Nothing},
    use::Nothing,
    options::Union{Dict{String,Any}, Nothing})
    return esc(
        quote
            Douglass.@generate_byvec!($t, $by, $sort, $arguments, $filter, $options )
        end
    )
end

# Stata's behavior:
#
#   - `if` restricts both the input and output rows to the ones that satisfy this condition
#   - `egen` does not by default restrict the input rows to those that are nonmissing for all RHS variables, but for some operations it does.
#       It never restricts the set of output rows to the nonmissing ones
#
# We follow Stata's behavor on [if], and, by default, restrict the input set to those that are nonmissing for all variables

# this is a new implementation that separates the parsing and variable generation from the parsing
macro generate_byvec!(t::Symbol,
    varlist_by::Union{Vector{Symbol}, Nothing},
    varlist_sort::Union{Vector{Symbol}, Nothing},
    arguments::Expr,
    filter::Union{Expr,Nothing},
    options::Union{Dict{String,Any}, Nothing})

    # assert that `arguments` is an assignment
    (arguments.head == :(.=)) && error("`egen` expects a vector-wise assignment operation, e.g. `:x = :y + :z`. Do not broadcast the assignment operator.")
    (arguments.head == :(=)) || error("`egen` expects an assignment operation, e.g. :x = :y + :z")

    assigned_var_type::Symbol = :Any
    if isexpr(arguments.args[1])
        # this needs to take the form :x::Type
        error("Expected a Symbol on the left-hand side of `egen` assignment operation, found Expr.")
        # ( arguments.args[1].head == Symbol("::") ) || error("Expected data type in assignment operation, e.g. `:x::Float64 = ...`")
        # (size(arguments.args[1].args,1) == 2 ) || error("Expected data type in assignment operation, e.g. `:x::Float64 = ...`")
        assigned_var_qn = arguments.args[1].args[1]::QuoteNode
    elseif isa(arguments.args[1], QuoteNode)
        assigned_var_qn = arguments.args[1]::QuoteNode
    end

    # if the RHS of the assignment expression is currently a symbol, make it an Expr
    #transformation = (typeof(arguments.args[2]) == Symbol) ? Expr(arguments.args[2]) : arguments.args[2]
    #if !isexpr(arguments.args[2])
    # transformation::Expr = :( )
    if isexpr(arguments.args[2])
        transformation = arguments.args[2]
    elseif isa(arguments.args[2], QuoteNode)
        transformation = Expr(:call, identity, arguments.args[2])
    else
        transformation = Expr(:call, identity, arguments.args[2])
    end

    # println("transformation is a $(typeof(transformation)) with value $(transformation).")
    # println("assigned_var_qn is a $(typeof(assigned_var_qn)) with value $(assigned_var_qn).")

    return esc(
        quote
            # check variable is not present
            ($(assigned_var_qn) ∉ propertynames($t)) || error("Variable $($(assigned_var_qn)) already present in DataFrame.")
            # do the assignment
            Douglass.@transform_byvec_allobs!($t, $varlist_by, $varlist_sort, $assigned_var_qn, $transformation, $filter)
        end
    )
end

# vector-wise transformation with `if` condition (`by` and `if`)
# reference implementation
macro transform_byvec_allobs!(t::Symbol,
    varlist_by::Vector{Symbol},
    varlist_sort::Union{Vector{Symbol}, Nothing},
    assigned_var_qn::QuoteNode,
    transformation,
    filter::Expr)

    # go through the expression tree and pick up QuoteNodes.
    qn_vec = Vector{QuoteNode}()
    transformation_converted = deepcopy(transformation)
    filter_converted = deepcopy(filter)

    # replace all QuoteNode's by args[x]
    replace_QuoteNodes!(transformation_converted, :args, qn_vec)
    # and the same for the filter condition
    replace_QuoteNodes!(filter_converted, :args, qn_vec)
    # add an index subscript to all arg's: arg[#] -> arg[#][index]
    # so that we can condition on a filter
    push_index!(transformation_converted, :indices, :args)

    # a vector of the Symbol's for each QuoteNode
    symbol_vec = [qn_vec[i].value for i=1:length(qn_vec)]

    # DEEBUGGING OUTPUT:
    # println("transformation_converted is a $(typeof(transformation_converted)) with value $(transformation_converted).")
    # println("filter_converted is a $(typeof(filter_converted)) with value $(filter_converted).")
    # println("qn_vec is a $(typeof(qn_vec)) with value $(qn_vec).")
    # println("first element is a $(typeof(qn_vec[1])) with value $(qn_vec[1]).")
    # @show qn_vec[1]
    # dump(transformation_converted)

    return esc(
        quote
            # sort, if we need to, (first by-variables, then sort-variables)
            if !isnothing($(varlist_sort)) && !isempty($(varlist_sort))
                sort!($t, vcat($varlist_by, $varlist_sort))
            end

            # define a barrier function, which will be compiled and type stable
            # see Ch. 11 of Bogumil Kaminski's DataFrames tutorial
            local function my_f(args...)
                # indices are the indices of the input columns to be used.
                indices = BitArray($filter_converted)
                # the following is a wrapper for a[indices] = $transformation_converted
                # but if $transformation_converted expands into a scalar, it's broadcast
                return Douglass.assign_helper_gen(indices,$transformation_converted)
            end
            transform!(groupby($t, $varlist_by, sort=false, skipmissing = true ), $(symbol_vec) => my_f => $(assigned_var_qn))
        end

    )
end

# `by` but no `if`
macro transform_byvec_allobs!(t::Symbol,
    varlist_by::Vector{Symbol},
    varlist_sort::Union{Vector{Symbol}, Nothing},
    assigned_var_qn::QuoteNode,
    transformation,
    filter::Nothing)

    # go through the expression tree and pick up QuoteNodes.
    qn_vec = Vector{QuoteNode}()
    transformation_converted = deepcopy(transformation)

    # replace all QuoteNode's by args[x]
    replace_QuoteNodes!(transformation_converted, :args, qn_vec)

    # a vector of the Symbol's for each QuoteNode
    symbol_vec = [qn_vec[i].value for i=1:length(qn_vec)]

    return esc(
        quote
            # sort, if we need to, (first by-variables, then sort-variables)
            if !isnothing($(varlist_sort)) && !isempty($(varlist_sort))
                sort!($t, vcat($varlist_by, $varlist_sort))
            end

            # define a barrier function, which will be compiled and type stable
            # see Ch. 11 of Bogumil Kaminski's DataFrames tutorial
            local function my_f(args...)
                return $transformation_converted
            end

            transform!(groupby($t, $varlist_by, sort=false, skipmissing = true ), $(symbol_vec) => my_f => $(assigned_var_qn))
        end

    )
end

# version without groups but with if
macro transform_byvec_allobs!(t::Symbol,
    varlist_by::Nothing,
    varlist_sort::Union{Vector{Symbol}, Nothing},
    assigned_var_qn::QuoteNode,
    transformation,
    filter::Expr)

    # go through the expression tree and pick up QuoteNodes.
    qn_vec = Vector{QuoteNode}()
    transformation_converted = deepcopy(transformation)
    filter_converted = deepcopy(filter)

    # replace all QuoteNode's by args[x]
    replace_QuoteNodes!(transformation_converted, :args, qn_vec)
    # and the same for the filter condition
    replace_QuoteNodes!(filter_converted, :args, qn_vec)
    # add an index subscript to all arg's: arg[#] -> arg[#][index]
    # so that we can condition on a filter
    push_index!(transformation_converted, :indices, :args)

    # a vector of the Symbol's for each QuoteNode
    symbol_vec = [qn_vec[i].value for i=1:length(qn_vec)]

    # println("transformation_converted is a $(typeof(transformation_converted)) with value $(transformation_converted).")
    # println("filter_converted is a $(typeof(filter_converted)) with value $(filter_converted).")
    # println("qn_vec is a $(typeof(qn_vec)) with value $(qn_vec).")
    # println("first element is a $(typeof(qn_vec[1])) with value $(qn_vec[1]).")
    # @show qn_vec[1]
    # dump(transformation_converted)

    return esc(
        quote
            # define a barrier function, which will be compiled and type stable
            # see Ch. 11 of Bogumil Kaminski's DataFrames tutorial
            local function my_f(args...)
                # indices are the indices of the input columns to be used.
                indices = BitArray($filter_converted)
                # the following is a wrapper for a[indices] = $transformation_converted
                # but if $transformation_converted expands into a scalar, it's broadcast
                return Douglass.assign_helper_gen(indices,$transformation_converted)
            end
            transform!($t, $(symbol_vec) => my_f => $(assigned_var_qn))
        end
    )
end
# version without `by` and `if`
macro transform_byvec_allobs!(t::Symbol,
    varlist_by::Nothing,
    varlist_sort::Union{Vector{Symbol}, Nothing},
    assigned_var_qn::QuoteNode,
    transformation,
    filter::Nothing)

    # go through the expression tree and pick up QuoteNodes.
    qn_vec = Vector{QuoteNode}()
    transformation_converted = deepcopy(transformation)

    # replace all QuoteNode's by args[x]
    replace_QuoteNodes!(transformation_converted, :args, qn_vec)

    # a vector of the Symbol's for each QuoteNode
    symbol_vec = [qn_vec[i].value for i=1:length(qn_vec)]

    return esc(
        quote
            # define a barrier function, which will be compiled and type stable
            # see Ch. 11 of Bogumil Kaminski's DataFrames tutorial
            local function my_f(args...)
                return $transformation_converted
            end
            transform!($t, $(symbol_vec) => my_f => $(assigned_var_qn))
        end
    )
end


text:
Base.@kwdef struct DGPMErlang <: AbstractModel
    # Data
    c0::Vector{Float64}
    y0::Vector{Float64}
    X0::Matrix{Float64}
    y1::Vector{Float64}
    X1::Matrix{Float64}
    # Transformed data
    event::Vector{Bool} = y0 .< c0
    N0::Int = length(y0)
    N1::Int = length(y1)
    # HyperParameters
    a_φ0::Float64 = 2.0
    b_φ0::Float64 = 0.1
    a_λ0::Float64 = 0.1
    b_λ0::Float64 = 0.1
    # Parameters
    φ::Vector{Float64} = rand(Gamma(a_φ0, 1.0 / b_φ0), 1)
    λ::Vector{Float64} = rand(Gamma(a_λ0, 1.0 / b_λ0), 1)
    # Transformed parameters
    sumy::Base.RefValue{Float64} = Ref(0.0)
    sumlogy::Vector{Float64} = zeros(1)
    S::Vector{Float64} = zeros(N1)
    # Skeleton
    skl::Skeleton = Skeleton(; y0, y1, X0, X1)
end

function skeleton(m::DGPMErlang)
    return m.skl
end

function kernel_pdf(m::DGPMErlang, yi::Float64, j::Int)
    kernel = Erlang(ceil(Int, m.φ[j]), 1.0 / m.λ[])
    return pdf(kernel, yi)
end

function kernel_cdf(m::DGPMErlang, yi::Float64, j::Int)
    kernel = Erlang(ceil(Int, m.φ[j]), 1.0 / m.λ[])
    return 1 - cdf(kernel, yi)
end

function update_atoms!(m::DGPMErlang)
    update_suffstats!(m)
    update_λ!(m)
    update_φ!(m)
    update_y!(m)
end

function update_suffstats!(m::DGPMErlang)
    (; y0, sumlogy, sumy, skl) = m
    d = cluster_labels(skl)
    if length(sumlogy) < rmax(skl)
        resize!(sumlogy, rmax(skl))
    end
    for j in 1:rmax(skl)
        sumlogy[j] = 0.0
    end
    sumy[] = 0.0
    for i in 1:length(y0)
        sumy[] += y0[i]
        sumlogy[d[i]] += log(y0[i])
    end
    return nothing
end

function update_λ!(m::DGPMErlang)
    (; λ, φ, a_λ0, b_λ0, sumy, skl) = m
    (; d) = skl
    a_λ1 = a_λ0
    b_λ1 = b_λ0 + sumy[]
     for di in d
        a_λ1 += ceil(φ[di])
    end
    λ[] = rand(Gamma(a_λ1, 1.0 / b_λ1))
    return nothing
end

function update_φ!(m::DGPMErlang)
    (; φ, a_φ0, b_φ0, skl) = m
    pφ = Gamma(a_φ0, 1.0 / b_φ0)
    while length(φ) < rmax(skl)
        push!(φ, rand(pφ))
    end
    for j in 1:rmax(skl)
        φ0 = φ[j]
        d0 = truncated(Normal(φ0, 0.1), 0, Inf)
        φ1 = rand(d0)
        d1 = truncated(Normal(φ1, 0.1), 0, Inf)
        log_ar = 0.0
        log_ar += (logpφ(m, φ1, j) + logpdf(d1, φ0))
        log_ar -= (logpφ(m, φ0, j) + logpdf(d0, φ1))
        if rand() < exp(log_ar)
            (φ[j] = φ1)
        end
    end
    return nothing
end

function logpφ(m::DGPMErlang, φ0, j)
    (; λ, a_φ0, b_φ0, sumlogy, skl) = m
    n = cluster_sizes(skl)
    δ0 = ceil(Int, φ0)
    return (
        n[j] * δ0 * log(λ[]) +
        (δ0 - 1) * sumlogy[j] -
        n[j] * logfactorial(δ0 - 1) +
        (a_φ0 - 1) * log(φ0) -
        b_φ0 * φ0
    )
end

function update_y!(m::DGPMErlang)
    (; y0, c0, event, φ, λ, skl) = m
    d = cluster_labels(skl)
    for i in 1:length(y0)
        if !event[i]
            pdf = Erlang(ceil(Int, φ[d[i]]), 1.0 / λ[])
            tpdf = Truncated(pdf, c0[i], Inf)
            if cdf(pdf, c0[i]) < 0.99
                y0[i] = rand(tpdf)
            end
        end
    end
end

function update_S!(m::DGPMErlang)
    (; S, skl) = m
    (; N1, y1, ϕ1, s, rmax) = skl
    for i in 1:N1
        S[i] = 0.0
        for j in 1:rmax[]
            wj = get_w(ϕ1[i], s[], j)
            S[i] += wj * kernel_cdf(m, y1[i], j)
        end
    end
end

function sample!(m::DGPMErlang; mcmcsize = 4000, burnin = 2000, thin = 1)
    (; S, skl) = m
    (; N1, D1, f, β) = skl
    chainf = [zeros(N1) for _ in 1:(mcmcsize - burnin) ÷ thin]
    chainS = [zeros(N1) for _ in 1:(mcmcsize - burnin) ÷ thin]
    chainβ = [zeros(D1) for _ in 1:(mcmcsize - burnin) ÷ thin]
    chaing = [zeros(Bool, D1) for _ in 1:(mcmcsize - burnin) ÷ thin]
    for iter in 1:mcmcsize
        step!(m)
        update_S!(m)
        if iter > burnin && iszero((iter - burnin) % thin)
            chainf[(iter - burnin) ÷ thin] .= f
            chainS[(iter - burnin) ÷ thin] .= S
            chainβ[(iter - burnin) ÷ thin] .= β
            chaing[(iter - burnin) ÷ thin] .= (β .!= 0.0)
        end
    end
    return chainf, chainS, chainβ, chaing
end


text:
module TempleMirror

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Temple Mirror" => Ahorn.EntityPlacement(
        Maple.TempleMirror,
        "rectangle"
    )
)

Ahorn.minimumSize(entity::Maple.TempleMirror) = 24, 24
Ahorn.resizable(entity::Maple.TempleMirror) = true, true

mirrorColor = (5, 7, 14, 1) ./ (255, 255, 255, 1)
mirrorTexture = "scenery/templemirror"

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

    width = Int(get(entity.data, "width", 48))
    height = Int(get(entity.data, "height", 32))

    return Ahorn.Rectangle(x, y, width, height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.TempleMirror, room::Maple.Room)
    width = Int(get(entity.data, "width", 48))
    height = Int(get(entity.data, "height", 32))

    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    frame = Ahorn.getSprite(mirrorTexture, "Gameplay")

    Ahorn.drawRectangle(ctx, 2, 2, width - 4, height - 4, mirrorColor)

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, frame, (i - 1) * 8, 0, 8, 0, 8, 8)
        Ahorn.drawImage(ctx, frame, (i - 1) * 8, height - 8, 8, 16, 8, 8)
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, 0, (i - 1) * 8, 0, 8, 8, 8)
        Ahorn.drawImage(ctx, frame, width - 8, (i - 1) * 8, 16, 8, 8, 8)
    end

    Ahorn.drawImage(ctx, frame, 0, 0, 0, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, width - 8, 0, 16, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, 0, height - 8, 0, 16, 8, 8)
    Ahorn.drawImage(ctx, frame, width - 8, height - 8, 16, 16, 8, 8)
end

end