text:
using MKL_jll, MKL, Test, SpecialFunctions
import LinearAlgebra

@show LinearAlgebra.BLAS.get_config()

@testset "Sanity Tests" begin
  @test LinearAlgebra.BLAS.get_config().loaded_libs[1].libname == libmkl_rt
  @test LinearAlgebra.peakflops() > 0
end

# Test #98 - issues with multi-threaded MKL on mac
# Test https://github.com/JuliaLang/julia/issues/40787

@testset "Threading nondeterminism test" begin
    A = [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -5 -5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -5 -5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -5 -5 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -5 -5 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -5 -5 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -5 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 -5 0 0 0 0 0 -5 -5 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 -5 0 0 0 0 0 -5 -5 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 -5 0 0 0 0 0 -5 -5 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 -5 0 0 0 0 0 -5 -5 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 -5 0 0 0 0 0 -5 -5 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 -5 0 0 0 0 0 -5 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 -5 0 0 0 0 0 -5 -5 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 -5 0 0 0 0 0 -5 -5 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 -5 0 0 0 0 0 -5 -5 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 -5 0 0 0 0 0 -5 -5 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 -5 0 0 0 0 0 -5 -5; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 -5 0 0 0 0 0 -5; -5 0 0 0 0 0 -5 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; -5 -5 0 0 0 0 0 -5 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 -5 -5 0 0 0 0 0 -5 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 -5 -5 0 0 0 0 0 -5 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 -5 -5 0 0 0 0 0 -5 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 -5 -5 0 0 0 0 0 -5 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 -5 0 0 0 0 0 -5 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 -5 -5 0 0 0 0 0 -5 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 -5 -5 0 0 0 0 0 -5 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 -5 -5 0 0 0 0 0 -5 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 -5 -5 0 0 0 0 0 -5 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 -5 -5 0 0 0 0 0 -5 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 -5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 -5 -5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 -5 -5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -5 -5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -5 -5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -5 -5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1];

    B = [1 0 0 0 0; 0 1 0 0 0; 0 0 1 0 0; 0 0 0 1 0; 0 0 0 0 1; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0];

    for i=1:10
        @test sum(abs, A \ B) ≈ 24.772054506344045
    end
end

# Run all the LinearAlgebra stdlib tests
include(joinpath(Sys.STDLIB, "LinearAlgebra", "test", "runtests.jl"))


text:
#!/usr/bin/env julia

function main(outlen:: Integer)
	chars = Vector{Char}(undef, outlen)
	for i in 1:outlen
		chars[i] = rand(Char)
	end
	print(String(chars))
end

main(parse(Int, ARGS[1]))

ext: jl

max_stars_repo_name: MasFlam/shell-utils

max_issues_repo_name: MasFlam/shell-utils

max_forks_repo_name: MasFlam/shell-utils


text:
"""

"""
function plotdset2Dsym(dset::ComsolUtils.CDataDSet;     # contains the dset data
        n=nothing,                                      # name of variable to be plotted
        t=nothing,                                      # time at which solution is plotted
        limits=nothing,                                 # limits [rmin rmax zmin zmax]
        datalim=nothing,                                # data limits [vmin vmax]
        levels=10,                                      # levels of contourf plot
        resolution=nothing,                             # resolution of figure
        N=100,                                          # number of grid points per dimension
        cmap=ColorSchemes.vik)

    ## PARSE INPUTS ========================================================= ##
    # parse name and unit
    n = isnothing(n) ? dset.n[1] : n
    unit = []
    for i in 1:length(dset.n); if n==dset.n[i]; unit=dset.u[i]; end; end

    # parse t
    It = isnothing(t) ? length(dset.t) : findfirst(x->x>t, dset.t)
    if isnothing(It); error("Cannot extrapolate time @ $t."); end

    # parse limits
    limits = isnothing(limits) ? [min(dset.p.r...), max(dset.p.r...), min(dset.p.z...), max(dset.p.z...)] : limits

    # extract data vector
    v = dset.v[n][:,It] |> Vector{Float64}

    # parse data limits
    datalim = isnothing(datalim) ? [min(filter(!isnan, v)...), max(filter(!isnan, v)...)] : datalim

    ## PREPARE DATA ========================================================= ##
    rgrid = LinRange(limits[1], limits[2], N) # grid vector
    zgrid = LinRange(limits[3], limits[4], N) # -
    rq = repeat(rgrid, N)       # query vector
    zq = repeat(zgrid, inner=N) # -

    ## INTERPOLATE & PLOT =================================================== ##
    vq = interpolation2Dsym(dset.p, v, rq, zq)
    V = reshape(vq, N, N)

    # create figure
    fig = isnothing(resolution) ? Figure() : Figure(resolution=resolution)
    ax = fig[1, 1] = Axis(fig)
    cf = contourf!(ax, rgrid, zgrid, V;
        levels=LinRange(datalim[1], datalim[2], levels),
        colormap=cmap,
        extendlow=cmap[1],
        extendhigh=cmap[end])
    Colorbar(fig[1,2], cf, tellheight=false, label="$(n) [$unit]", ticks=LinRange(datalim[1], datalim[2], levels))

    # format axis & add labels
    limits!(ax, limits...)
    ax.aspect = DataAspect();
    ax.xlabel = "r coordinate [m]"
    ax.ylabel = "z coordinate [m]"

    # return figure handle
    return fig
end

## INTERNAL FUNCTIONS ======================================================= ##
function nicelabel(label::String)
    substitutions = [
        ## superscripts
        "^0"=>"⁰", "^1"=>"¹", "^2"=>"²", "^3"=>"³", "^4"=>"⁴", "^5"=>"⁵",
        "^6"=>"⁶", "^7"=>"⁷", "^8"=>"⁸", "^9"=>"⁹",

        # q - missing
        "^a"=>"ᵃ", "^b"=>"ᵇ", "^c"=>"ᶜ", "^d"=>"ᵈ", "^e"=>"ᵉ", "^f"=>"ᶠ",
        "^g"=>"ᵍ", "^h"=>"ʰ", "^i"=>"ⁱ", "^j"=>"ʲ", "^k"=>"ᵏ", "^l"=>"ˡ",
        "^m"=>"ᵐ", "^n"=>"ⁿ", "^o"=>"ᵒ", "^p"=>"ᵖ", "^r"=>"ʳ", "^s"=>"ˢ",
        "^t"=>"ᵗ", "^u"=>"ᵘ", "^v"=>"ᵛ", "^w"=>"ʷ", "^x"=>"ˣ", "^y"=>"ʸ",
        "^z"=>"ᶻ",

        ## subscripts
        "_0"=>"₀", "_1"=>"₁", "_2"=>"₂", "_3"=>"₃", "_4"=>"₄", "_5"=>"₅",
        "_6"=>"₆", "_7"=>"₇", "_8"=>"₈", "_9"=>"₉",

        # b, c, d, f, g, q, w, y, z - missing
        "_a"=>"ₐ", "_e"=>"ₑ", "_h"=>"ₕ", "_i"=>"ᵢ", "_j"=>"ⱼ", "_k"=>"ₖ",
        "_l"=>"ₗ", "_m"=>"ₘ", "_n"=>"ₙ", "_o"=>"ₒ", "_p"=>"ₚ", "_r"=>"ᵣ",
        "_s"=>"ₛ", "_t"=>"ₜ",  "_u"=>"ᵤ", "_v"=>"ᵥ", "_x"=>"ₓ"
        ]

    for sub in substitutions
        label = replace(label, sub)
    end

    return label
end


text:
#=
Authors:
Darius Russell Kish

Created: April 14, 2020
Modified: April 14, 2020

This implementation is based strongly in the Disjoint Set implementation
in JuliaCollections' DataStructures which is under active development. It
represents proper implementation of data structures for Julia and much of the
typing and integration into the language is taken from their implementation,
however the function implementations are my own.
https://github.com/JuliaCollections/DataStructures.jl/blob/master/src/disjoint_set.jl

Implementations for collections of arbitrary inputs are *not* provided here
and it is assumed this low-level data structure will be wrapped for such
functionality. This avoids issues of key-value association time complexity for
conversion between elements and their index value in analyses of this data
structure.
=#

#=
Notation T <: Integer :: T is a subtype of Integer

An example type tree:
abstract type Number end
abstract type Real     <: Number end
abstract type AbstractFloat <: Real end
abstract type Integer  <: Real end
abstract type Signed   <: Integer end
abstract type Unsigned <: Integer end
=#
mutable struct QuickFind{T<:Integer}
    ids::Vector{T}
    nGroups::T
end

# Construction type 1
# Notation: collect(Base.OneTo(n)) == list(range(n)) in Python
QuickFind(n::T) where {T <: Integer} =
    QuickFind{T}(collect(Base.OneTo(n)), n)

# Construction type 2
QuickFind{T}(n::Integer) where {T <: Integer} =
    QuickFind{T}(collect(Base.OneTo(n)), n)

# The length or "size" of a UnionFind is the length of either of its vectors
length(uf::QuickFind) = Base.length(uf.ids)

"""
    num_groups(uf::QuickFind)
Get a number of groups.
"""
num_groups(uf::QuickFind) = uf.nGroups

"""
    find(uf::QuickFind{T}, x::T, y::T)
Determine if x and y belong to the same group
"""
function find(uf::QuickFind{T}, x::T, y::T) where {T <: Integer}
    return uf.ids[x] == uf.ids[y]
end

"""
    find(uf::QuickFind{T}, t::Tuple{T,T})
Determine if x and y belong to the same group
"""
function find(uf::QuickFind{T}, t::Tuple{T,T}) where {T <: Integer}
    return uf.ids[t[1]] == uf.ids[t[2]]
end

"""
    unite(uf::QuickFind{T}, x::T, y::T)
Merges two groups together
"""
function unite!(uf::QuickFind{T}, x::T, y::T) where {T <: Integer}
    # We cannot guarantee x,y are in bounds
    pid = uf.ids[x]
    tid = uf.ids[y]
    # We can guarantee the below loop is in bounds, so we
    # can turn bounds checking off
    @inbounds for i = 1:length(uf)
        if uf.ids[i] == pid
            uf.ids[i] = tid
        end
    end
    # Groups are merged so decrement
    uf.nGroups -= 1
end

"""
    unite(uf::QuickFind{T}, t::Tuple{T,T})
Merges two groups together
"""
function unite!(uf::QuickFind{T}, t::Tuple{T,T}) where {T <: Integer}
    # We cannot guarantee x,y are in bounds
    pid = uf.ids[t[1]]
    tid = uf.ids[t[2]]
    # We can guarantee the below loop is in bounds, so we
    # can turn bounds checking off
    @inbounds for i = 1:length(uf)
        if uf.ids[i] == pid
            uf.ids[i] = tid
        end
    end
    # Groups are merged so decrement
    uf.nGroups -= 1
end


text:
using Distances
using NearestNeighbors
using StaticArrays

export NeighborList, neighborlist

struct NeighborList
    j::Vector{Vector{Int}}
    R::Vector{Vector{Float64}}
    r::Vector{Vector{SVector{3,Float64}}}
end
Base.length(nn::NeighborList) = length(nn.j)

function get_distance(L::SVector{3,<:AbstractFloat}, x::SVector{3,<:AbstractFloat}, y::SVector{3,<:AbstractFloat})
    broadcast(L, x, y) do Li, xi, yi
        if Li == Inf
            xi - yi
        else
            d = mod(xi - yi, Li)
            2d < Li ? d : d - Li
        end
    end
end

function neighborlist(A::AbstractSystem{3}, rcutoff::Float64)
    # Convert Positions to Matrix for Ball tree
    X = [SVector{3}(ustrip.(p)) for p ∈ position(A)]

    # Create Metric for Periodic Boundary Conditions
    periodic = periodicity(A)
    L = @SVector [periodic[i] ? ustrip(bounding_box(A)[i][i]) : Inf for i = 1:3]
    d = Distances.PeriodicEuclidean(L)

    # Build Ball tree
    tree = BallTree(X, d)

    # Intialize empty vectors
    j = Vector{Vector{Int64}}(undef, length(X))              # Neighbors
    R = Vector{Vector{Float64}}(undef, length(X))            # Distances
    r = Vector{Vector{SVector{3,Float64}}}(undef, length(X)) # Positions

    # Fill vectors
    for n in 1:length(X)
        neighbors = filter(m -> m > n, inrange(tree, X[n], rcutoff, true))
        j[n] = zeros(Int64, length(neighbors))
        R[n] = zeros(Float64, length(neighbors))
        r[n] = zeros(SVector{3,Float64}, length(neighbors))
        for (i, m) in enumerate(neighbors)
            rr = get_distance(L, X[n], X[m])
            j[n][i] = m
            R[n][i] = norm(rr)
            r[n][i] = rr
        end
    end
    NeighborList(j, R, r)
end


text:
@testset "Gauss–Radau" begin
    n = 1
    x, w = gaussradau(n)
    @test x[1] ≈ -1.
    @test w[1] ≈ 2.

    n = 2
    x, w = gaussradau(n)
    @test x[1] ≈ -1.
    @test x[2] ≈ 1/3
    @test w[1] ≈ .5
    @test w[2] ≈ 1.5

    tol = 1e-14
    n = 42
    ntests = 4
    x, w = gaussradau(n)
    @test length(x) == n && length(w) == n
    @test abs(x[37] - 0.908847278001044) < tol
    @test abs(w[37] - 0.031190846817016) < tol
    @test x[1] == -1

    @test dot( w,exp.(x)) ≈ exp(1)-exp(-1)
end


text:
using ChannelAttribution
using Test
using Random

@test isempty(detect_ambiguities(ChannelAttribution))

include("heuristic.jl")


text:
using MetacommunityDynamics
using DynamicGrids
using Plots
using Distributions

# Spatially explicit levins model

rs = SpatiallyExplicitLevinsColonization{:O}(probability=0.3) +
     AbioticExtinction{:O}(:A, baseprobability=0.1);

pops = pointstogrid(generate(PoissonProcess, 20), gridsize=100)

areas = generate(StaticEnvironmentalLayer, similar(pops), Exponential(3), mask=pops)

areas = rand(5,5)

initocc = rand(Occupancy, 0.8, 5,5)
arrayout = ArrayOutput((O=initocc, A=areas ), tspan=1:3)
sim!(arrayout, rs)


text:
"""
$(TYPEDEF)

Epigraph of real or complex infinity norm cone of dimension `dim`.

    $(FUNCTIONNAME){T, R}(dim::Int, use_dual::Bool = false)
"""
mutable struct EpiNormInf{T <: Real, R <: RealOrComplex{T}} <: Cone{T}
    use_dual_barrier::Bool
    dim::Int
    d::Int
    is_complex::Bool

    point::Vector{T}
    dual_point::Vector{T}
    grad::Vector{T}
    dder3::Vector{T}
    vec1::Vector{T}
    vec2::Vector{T}
    feas_updated::Bool
    grad_updated::Bool
    hess_updated::Bool
    inv_hess_updated::Bool
    inv_hess_aux_updated::Bool
    is_feas::Bool
    hess::Symmetric{T, SparseMatrixCSC{T, Int}}
    inv_hess::Symmetric{T, Matrix{T}}

    w::Vector{R}
    mu::Vector{R}
    zeta::Vector{T}
    cu::T
    Zu::T
    wumzi::Vector{R}
    zti::Vector{T}

    s1::Vector{T}
    s2::Vector{T}
    w1::Vector{R}
    w2::Vector{R}

    function EpiNormInf{T, R}(
        dim::Int;
        use_dual::Bool = false,
        ) where {T <: Real, R <: RealOrComplex{T}}
        @assert dim >= 2
        cone = new{T, R}()
        cone.use_dual_barrier = use_dual
        cone.dim = dim
        cone.is_complex = (R <: Complex)
        @assert !cone.is_complex || iseven(dim - 1)
        cone.d = (cone.is_complex ? div(dim - 1, 2) : dim - 1)
        return cone
    end
end

reset_data(cone::EpiNormInf) = (cone.feas_updated = cone.grad_updated =
    cone.hess_updated = cone.inv_hess_updated = cone.inv_hess_aux_updated = false)

use_sqrt_hess_oracles(::Int, cone::EpiNormInf) = false

function setup_extra_data!(
    cone::EpiNormInf{T, R},
    ) where {T <: Real, R <: RealOrComplex{T}}
    d = cone.d
    cone.w = zeros(R, d)
    cone.mu = zeros(R, d)
    cone.zeta = zeros(T, d)
    cone.wumzi = zeros(R, d)
    cone.s1 = zeros(T, d)
    cone.s2 = zeros(T, d)
    cone.w1 = zeros(R, d)
    if cone.is_complex
        cone.w2 = zeros(R, d)
    else
        cone.zti = zeros(T, d)
    end
    return cone
end

get_nu(cone::EpiNormInf) = 1 + cone.d

function set_initial_point!(
    arr::AbstractVector{T},
    cone::EpiNormInf{T},
    ) where {T <: Real}
    arr .= 0
    arr[1] = sqrt(T(get_nu(cone)))
    return arr
end

function update_feas(cone::EpiNormInf{T}) where T
    @assert !cone.feas_updated
    u = cone.point[1]

    if u > eps(T)
        @views vec_copyto!(cone.w, cone.point[2:end])
        cone.is_feas = (u - maximum(abs, cone.w) > eps(T))
    else
        cone.is_feas = false
    end

    cone.feas_updated = true
    return cone.is_feas
end

function is_dual_feas(cone::EpiNormInf{T}) where T
    dp = cone.dual_point
    u = dp[1]
    if u > eps(T)
        if cone.is_complex
            @inbounds norm1 = sum(hypot(dp[2i], dp[2i + 1]) for i in 1:cone.d)
        else
            @views norm1 = sum(abs, dp[2:end])
        end
        return (u - norm1 > eps(T))
    end
    return false
end

remul(a::T, b::T) where {T <: Real} = a * b

remul(a::Complex{T}, b::Complex{T}) where {T <: Real} =
    (real(a) * real(b) + imag(a) * imag(b))

function update_grad(cone::EpiNormInf{T}) where T
    @assert cone.is_feas
    u = cone.point[1]
    w = cone.w
    mu = cone.mu
    zeta = cone.zeta
    g = cone.grad

    @. mu = w / u
    @. zeta = T(0.5) * (u - remul(mu, w))
    cone.cu = (cone.d - 1) / u

    g[1] = cone.cu - sum(inv, zeta)

    if cone.is_complex
        @. cone.w1 = mu / zeta
        @views vec_copyto!(g[2:end], cone.w1)
    else
        @. g[2:end] = mu / zeta
    end

    cone.grad_updated = true
    return cone.grad
end

function update_hess(cone::EpiNormInf)
    isdefined(cone, :hess) || alloc_hess!(cone)
    d = cone.d
    u = cone.point[1]
    zeta = cone.zeta
    tzi2 = cone.s1
    g = cone.grad
    Hnz = cone.hess.data.nzval # modify nonzeros of upper triangle

    ui = inv(u)
    @. tzi2 = (inv(zeta) - ui) / zeta
    Hnz[1] = sum(tzi2) - cone.cu / u

    k = 2
    if cone.is_complex
        @inbounds for i in 1:d
            z_i = zeta[i]
            (mur, mui) = reim(cone.mu[i])
            mzir = g[2i]
            mzii = g[2i + 1]
            Hnz[k] = -mzir / z_i
            Hnz[k + 1] = (mur * mzir + ui) / z_i
            Hnz[k + 2] = -mzii / z_i
            Hnz[k + 3] = mzir * mzii
            Hnz[k + 4] = (mui * mzii + ui) / z_i
            k += 5
        end
    else
        @inbounds for i in 1:d
            Hnz[k] = -g[1 + i] / zeta[i]
            Hnz[k + 1] = tzi2[i]
            k += 2
        end
    end

    cone.hess_updated = true
    return cone.hess
end

function hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::EpiNormInf{T, T},
    ) where T
    @assert cone.grad_updated
    u = cone.point[1]
    mu = cone.mu
    zeta = cone.zeta
    s1 = cone.s1

    @inbounds for j in 1:size(prod, 2)
        p = arr[1, j]
        @views r = arr[2:end, j]

        pui = p / u
        @. s1 = (p - mu * r) / zeta

        prod[1, j] = sum((s1[i] - pui) / zeta[i] for i in 1:cone.d) -
            cone.cu * pui

        @. prod[2:end, j] = (r / u - s1 * mu) / zeta
    end

    return prod
end

function hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::EpiNormInf{T, Complex{T}},
    ) where T
    @assert cone.grad_updated
    u = cone.point[1]
    mu = cone.mu
    zeta = cone.zeta
    r = cone.w1
    w2 = cone.w2
    s1 = cone.s1

    @inbounds for j in 1:size(prod, 2)
        p = arr[1, j]
        @views vec_copyto!(r, arr[2:end, j])

        pui = p / u
        @. s1 = (p - remul(mu, r)) / zeta

        prod[1, j] = sum((s1[i] - pui) / zeta[i] for i in 1:cone.d) -
            cone.cu * pui

        @. w2 = (r / u - s1 * mu) / zeta
        @views vec_copyto!(prod[2:end, j], w2)
    end

    return prod
end

function update_inv_hess_aux(cone::EpiNormInf)
    @assert !cone.inv_hess_aux_updated
    @assert cone.grad_updated
    u = cone.point[1]
    w = cone.w
    wumzi = cone.wumzi
    umz = cone.s1

    @. umz = u - cone.zeta
    cone.Zu = -cone.cu + sum(inv, umz)

    @. wumzi = w / umz

    if !cone.is_complex
        @. cone.zti = u - wumzi * w
    end

    cone.inv_hess_aux_updated = true
end

function update_inv_hess(cone::EpiNormInf)
    cone.inv_hess_aux_updated || update_inv_hess_aux(cone)
    isdefined(cone, :inv_hess) || alloc_inv_hess!(cone)
    d = cone.d
    u = cone.point[1]
    w = cone.w
    zeta = cone.zeta
    wumzi = cone.wumzi
    Hi = cone.inv_hess.data

    hiuu = Hi[1, 1] = u / cone.Zu

    @views Hiuw = Hi[1, 2:end]
    @views Hiuw2 = Hi[2:end, 1]
    vec_copyto!(Hiuw2, wumzi)
    @. Hiuw = Hiuw2 * hiuu
    @views mul!(Hi[2:end, 2:end], Hiuw2, Hiuw')

    if cone.is_complex
        @inbounds for i in 1:d
            z_i = zeta[i]
            (wr, wi) = reim(w[i])
            (wumzir, wumzii) = reim(wumzi[i])
            k = 2 * i
            k1 = k + 1
            Hi[k, k] += z_i * (u - wumzir * wr)
            Hi[k1, k1] += z_i * (u - wumzii * wi)
            Hi[k, k1] -= z_i * wumzir * wi
        end
    else
        @inbounds for i in 1:d
            k = 1 + i
            Hi[k, k] += zeta[i] * cone.zti[i]
        end
    end

    cone.inv_hess_updated = true
    return cone.inv_hess
end

function inv_hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::EpiNormInf{T, T},
    ) where T
    cone.inv_hess_aux_updated || update_inv_hess_aux(cone)
    u = cone.point[1]
    wumzi = cone.wumzi

    @inbounds for j in 1:size(prod, 2)
        p = arr[1, j]
        @views r = arr[2:end, j]

        c1 = u * (p + dot(wumzi, r)) / cone.Zu
        prod[1, j] = c1

        @. prod[2:end, j] = c1 * wumzi + cone.zeta * r * cone.zti
    end

    return prod
end

function inv_hess_prod!(
    prod::AbstractVecOrMat,
    arr::AbstractVecOrMat,
    cone::EpiNormInf{T, Complex{T}},
    ) where T
    cone.inv_hess_aux_updated || update_inv_hess_aux(cone)
    u = cone.point[1]
    w = cone.w
    wumzi = cone.wumzi
    s1 = cone.s1
    r = cone.w1
    w2 = cone.w2

    @inbounds for j in 1:size(prod, 2)
        p = arr[1, j]
        @views vec_copyto!(r, arr[2:end, j])

        @. s1 = remul(wumzi, r)

        c1 = u * (p + sum(s1)) / cone.Zu
        prod[1, j] = c1

        @. w2 = c1 * wumzi + cone.zeta * (u * r - s1 * w)
        @views vec_copyto!(prod[2:end, j], w2)
    end

    return prod
end

function dder3(cone::EpiNormInf{T, T}, dir::AbstractVector{T}) where T
    @assert cone.grad_updated
    u = cone.point[1]
    mu = cone.mu
    zeta = cone.zeta
    rui = cone.w1
    s1 = cone.s1
    s2 = cone.s2
    dder3 = cone.dder3

    p = dir[1]
    @views r = dir[2:end]

    pui = p / u
    @. rui = r / u
    @. s1 = (p - mu * r) / zeta
    @. s2 = 0.5 * (p * pui - rui * r) / zeta - abs2(s1)

    @inbounds c1 = sum((s1[i] * pui + s2[i]) / zeta[i] for i in 1:cone.d)
    dder3[1] = -c1 - cone.cu * abs2(pui)

    @. dder3[2:end] = (s1 * rui + s2 * mu) / zeta

    return dder3
end

function dder3(cone::EpiNormInf{T, Complex{T}}, dir::AbstractVector{T}) where T
    @assert cone.grad_updated
    u = cone.point[1]
    mu = cone.mu
    zeta = cone.zeta
    rui = cone.w1
    r = w2 = cone.w2
    s1 = cone.s1
    s2 = cone.s2
    dder3 = cone.dder3

    p = dir[1]
    @views vec_copyto!(r, dir[2:end])

    pui = p / u
    @. rui = r / u
    @. s1 = (p - remul(mu, r)) / zeta
    @. s2 = 0.5 * (p * pui - remul(rui, r)) / zeta - abs2(s1)

    @inbounds c1 = sum((s1[i] * pui + s2[i]) / zeta[i] for i in 1:cone.d)
    dder3[1] = -c1 - cone.cu * abs2(pui)

    @. w2 = (s1 * rui + s2 * mu) / zeta
    @views vec_copyto!(dder3[2:end], w2)

    return dder3
end

function alloc_hess!(cone::EpiNormInf{T, T}) where {T <: Real}
    # initialize sparse idxs for upper triangle of Hessian
    dim = cone.dim
    nnz_tri = 2 * dim - 1
    I = Vector{Int}(undef, nnz_tri)
    J = Vector{Int}(undef, nnz_tri)
    idxs1 = 1:dim
    @views I[idxs1] .= 1
    @views J[idxs1] .= idxs1
    idxs2 = (dim + 1):(2 * dim - 1)
    @views I[idxs2] .= 2:dim
    @views J[idxs2] .= 2:dim
    V = ones(T, nnz_tri)
    cone.hess = Symmetric(sparse(I, J, V, dim, dim), :U)
    return
end

function alloc_hess!(cone::EpiNormInf{T, Complex{T}}) where {T <: Real}
    # initialize sparse idxs for upper triangle of Hessian
    dim = cone.dim
    nnz_tri = 2 * dim - 1 + cone.d
    I = Vector{Int}(undef, nnz_tri)
    J = Vector{Int}(undef, nnz_tri)
    idxs1 = 1:dim
    @views I[idxs1] .= 1
    @views J[idxs1] .= idxs1
    idxs2 = (dim + 1):(2 * dim - 1)
    @views I[idxs2] .= 2:dim
    @views J[idxs2] .= 2:dim
    idxs3 = (2 * dim):nnz_tri
    @views I[idxs3] .= 2:2:dim
    @views J[idxs3] .= 3:2:dim
    V = ones(T, nnz_tri)
    cone.hess = Symmetric(sparse(I, J, V, dim, dim), :U)
    return
end

# TODO remove this in favor of new hess_nz_count etc functions
# that directly use uu, uw, ww etc
hess_nz_count(cone::EpiNormInf{<:Real, <:Real}) =
    3 * cone.dim - 2

hess_nz_count(cone::EpiNormInf{<:Real, <:Complex}) =
    3 * cone.dim - 2 + 2 * cone.d

hess_nz_count_tril(cone::EpiNormInf{<:Real, <:Real}) =
    2 * cone.dim - 1

hess_nz_count_tril(cone::EpiNormInf{<:Real, <:Complex}) =
    2 * cone.dim - 1 + cone.d

hess_nz_idxs_col(cone::EpiNormInf{<:Real, <:Real}, j::Int) =
    (j == 1 ? (1:cone.dim) : [1, j])

hess_nz_idxs_col(cone::EpiNormInf{<:Real, <:Complex}, j::Int) =
    (j == 1 ? (1:cone.dim) : (iseven(j) ? [1, j, j + 1] : [1, j - 1, j]))

hess_nz_idxs_col_tril(cone::EpiNormInf{<:Real, <:Real}, j::Int) =
    (j == 1 ? (1:cone.dim) : [j])

hess_nz_idxs_col_tril(cone::EpiNormInf{<:Real, <:Complex}, j::Int) =
    (j == 1 ? (1:cone.dim) : (iseven(j) ? [j, j + 1] : [j]))


text:
using LinearAlgebra, LowRankApprox, Formatting, ImageFiltering, Images, Plots, ImageMagick, Colors, TestImages, ImageView, ImageTransformations

I=load("rgb.jpg");

G=Gray.(I);
# size(G)

A=imresize(G, ratio=1/3);
# rank(A), size(A), eltype(A)

(p,q)=size(A);
B=Array{Float64}(A);

using DSP, Wavelets
B1=dwt(B, wavelet(WT.haar), 1);
B2=dwt(B, wavelet(WT.haar), 2);
B3=dwt(B, wavelet(WT.haar), 3);
B4=dwt(B, wavelet(WT.haar), 4);

B1a=zeros(Float64, p, q);
B2a=zeros(Float64, p, q);
B3a=zeros(Float64, p, q);
B4a=zeros(Float64, p, q);

B1a[1:Int(p/2),1:Int(q/2)]=B1[1:Int(p/2),1:Int(q/2)];
B2a[1:Int(p/4),1:Int(q/4)]=B2[1:Int(p/4),1:Int(q/4)];
B3a[1:Int(p/8),1:Int(q/8)]=B3[1:Int(p/8),1:Int(q/8)];
B4a[1:Int(p/16),1:Int(q/16)]=B4[1:Int(p/16),1:Int(q/16)];
B1d=B1-B1a;
B2d=B2-B2a;
B3d=B3-B3a;
B4d=B4-B4a;

C1a=idwt(B1a, wavelet(WT.haar), 1);
C2a=idwt(B2a, wavelet(WT.haar), 2);
C3a=idwt(B3a, wavelet(WT.haar), 3);
C4a=idwt(B4a, wavelet(WT.haar), 4);
C1d=idwt(B1d, wavelet(WT.haar), 1);
C2d=idwt(B2d, wavelet(WT.haar), 2);
C3d=idwt(B3d, wavelet(WT.haar), 3);
C4d=idwt(B4d, wavelet(WT.haar), 4);

A1a=Gray.(C1a);
A1d=Gray.(C1d);
A2a=Gray.(C2a);
A2d=Gray.(C2d);
A3a=Gray.(C3a);
A3d=Gray.(C3d);
A4a=Gray.(C4a);
A4d=Gray.(C4d);

plot(Gray.(A),
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="Level 0: Original Grayscale Image")
savefig("grayscale_haar_level0.png")

plot(Gray.(B1),
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="Level 1 (haar)")
savefig("grayscale_haar_level1.png")

plot(Gray.(C1a),
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="Level 1 Approximation (haar)")
savefig("grayscale_haar_level1_a.png")

plot(Gray.(C2a),
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="Level 2 Approximation (haar)")
savefig("grayscale_haar_level2_a.png")

plot(Gray.(C3a),
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="Level 3 Approximation (haar)")
savefig("grayscale_haar_level3_a.png")

plot(Gray.(C4a),
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="Level 4 Approximation (haar)")
savefig("grayscale_haar_level4_a.png")

plot(Gray.(C4d),
    xaxis=false,
    xticks=false,
    yaxis=false,
    yticks=false,
    grid=false,
    title="Level 4 Detail (haar)")
savefig("grayscale_haar_level4_d.png")


text:
# this is the 2nd attempt; in this attempt I will try use MCTS and gradients
using StatsBase, DataFrames, Flux, JLD2, Plots, JLD2, FileIO, BSON
include("2048_fn.jl")

policy_loss(x,y) = Flux.crossentropy(policy(x), y)
critic_loss(x,y) = sum((critic(x).-y).^2)

policy = Flux.Chain(
  x -> reshape(Float64.(x), (4,4,1,1)),
  Conv((2,2), 1=>64, relu),
  x -> meanpool(x, (2,2)),
  x -> reshape(x, (64,)),
  Dense(64, 32, relu),
  Dense(32, 4, relu),
  softmax);

critic = Flux.Chain(
  gridmove -> (reshape(Float64.(gridmove[1]), (4,4,1,1)), Float64.(Flux.onehot(gridmove[2], DIRS))),
  gridmove -> (Conv((2,2),1=>64, relu)(gridmove[1]) , gridmove[2]),
  gridmove -> (meanpool(gridmove[1], (2,2)), gridmove[2]),
  gridmove -> vcat(reshape(gridmove[1], (64,)), gridmove[2]),
  Dense(64+4, 32, relu),
  Dense(32, 1, relu)
  );
#
# policy_weights = Tracker.data.(params(policy));
# critic_weights = Tracker.data.(params(critic));
#
# BSON.@save "actor_critic_model1.bson" policy
# BSON.@load "actor_critic_model1.bson" policy
# BSON.@save "actor_critic_model2.bson" policy
#
# BSON.@load "actor_critic_model.bson" policy_weights critic_weights
# Flux.loadparams!(policy, policy_weights)
# Flux.loadparams!(critic, critic_weights)

opt = ADAM(params(policy))
opt_critic = ADAM(params(critic))

function train1!(policy, critic, policy_loss, critic_loss, opt, opt_critic, n;
  γ = 0.99 # next critic moves downward rate
  )
  for i = 1:n
    agrid = init_game()
    vm = valid_moves(agrid)
    while length(vm) > 0
      # use policy to generate searches
      p = Tracker.data(policy(agrid))[indexin(vm, DIRS)]

      sampled_move = sample(vm, Weights(p))
      (new_grid, _, _, reward) = simulate_one!(copy(agrid), sampled_move)

      # sample a new move for the new_grid
      new_vm = valid_moves(new_grid)

      if length(new_vm) > 0
        p = Tracker.data(policy(new_grid))[indexin(new_vm, DIRS)]
        sampled_next_action = sample(new_vm, Weights(p))

        # evaluate critic at the new location
        discounted_critic_value = γ*Tracker.data(critic((new_grid, sampled_next_action)))[1]
      else
        discounted_critic_value = 0
      end

      # update the critic
      Flux.train!(critic_loss, [((agrid, sampled_move), reward + discounted_critic_value)], opt_critic)

      # now use the updated critic to update to the policy
      best_move_according_to_new_critic = findmax(Dict{Symbol, Float64}((vm1, Tracker.data(critic((agrid, vm1)))[1]) for vm1 in vm))[2]
      Flux.train!(policy_loss, [(agrid, Flux.onehot(best_move_according_to_new_critic, DIRS))], opt)

      agrid = new_grid
      vm = new_vm
    end
    i += 1
    if i % 1000 == 0
      @save "actor_critic_model.jld2" policy critic
      display(agrid)
    end
  end
end

init = init_game()

policy(init)
[critic((init,vm)) for vm in DIRS]
@time train1!(policy, critic, policy_loss, critic_loss, opt, opt_critic, 1_000_000)
policy(init)
critic((init,:left))

(init, fnl, seq, cart, tf, rewards) = simulate_game!(init)

nmovesback = 1
last_grid = sim_seq(init, seq[1:end-nmovesback], cart[1:end-nmovesback], tf[1:end-nmovesback])
display(last_grid)
display(policy(last_grid))
display([(vm, Tracker.data(critic((last_grid, vm)))[1]) for vm in valid_moves(last_grid)])

#manually train this
agrid = copy(last_grid)
vm = valid_moves(agrid)

# use policy to generate searches
p = Tracker.data(policy(agrid))[indexin(vm, DIRS)]

sampled_move = sample(vm, Weights(p))
(new_grid, _, _, reward) = simulate_one!(copy(agrid), sampled_move)

# sample a new move for the new_grid
new_vm = valid_moves(new_grid)

if length(new_vm) > 0
  p = Tracker.data(policy(new_grid))[indexin(new_vm, DIRS)]
  sampled_next_action = sample(new_vm, Weights(p))

  # evaluate critic at the new location
  discounted_critic_value = γ*Tracker.data(critic((new_grid, sampled_next_action)))[1]
else
  discounted_critic_value = 0
end

critic((agrid, :right))
# update the critic
Flux.train!(critic_loss, [((agrid, sampled_move), reward + discounted_critic_value)], opt_critic)
display((critic((agrid, :right))|> Tracker.data, critic((agrid, :left))|>Tracker.data))

# now use the updated critic to update to the policy
best_move_according_to_new_critic = findmax(Dict{Symbol, Float64}((vm1, Tracker.data(critic((agrid, vm1)))[1]) for vm1 in vm))[2]
Flux.train!(policy_loss, [(agrid, Flux.onehot(best_move_according_to_new_critic, DIRS))], opt)

agrid = new_grid
old_grid = copy(new_grid)
display(agrid)
vm = new_vm

critic((init,:left))

agrid = last_grid

vm = valid_moves(agrid)

(init, fnl, seq, cart, tf, reward) = simulate_game!(agrid)

# obtain the last move
last_move = 2
before_last = sim_seq(init, seq[1:end-last_move], cart[1:end-last_move], tf[1:end-last_move])

vm = valid_moves(before_last)
@time avg_rewards_next_move(before_last, 100)

function train_n_games(n, model, opt)
  for i=1:n
    agrid = init_game()

    vm = valid_moves(agrid)
    while length(vm) > 0
      # compute average rewards
      avgr = avg_rewards_next_move(agrid, 100)
      grid_as_features = grid_to_feature(agrid)
      #model(grid_as_features)
      Flux.train!(loss, [(grid_as_features, findmax(avgr)[2])], opt)

      # now choose an action based on the updated policy
      themove = sample(vm, Weights(Tracker.data(model(grid_as_features))[indexin([:left,:up], DIRS)]))
      simulate_one!(agrid, themove)
      vm = valid_moves(agrid)
    end
    println(i)
    display(agrid)

    @save "policy_gradient_model.jld2" model

  end
end

@time train_n_games(1, model, opt)

agrid = init_game()
(init,_,seq,cart, tf,_) = simulate_game!(agrid)

new_grid = sim_seq(init, seq[1:end-2], cart[1:end-2], tf[1:end-2])

model(grid_to_feature(new_grid))
avgr = avg_rewards_next_move(new_grid, 100)
best_move = findmax(avgr)[2]
Flux.@epochs 10000 Flux.train!(loss, [(grid_to_feature(new_grid), Flux.onehot(best_move, DIRS))], opt)
model(grid_to_feature(new_grid))

Float64.(Flux.onehot(best_move, DIRS))

crossentropy(model(grid_to_feature(new_grid)), Float64.(Flux.onehot(best_move, DIRS)))

logloss(model(grid_to_feature(new_grid)), Float64.(Flux.onehot(:up, DIRS)))

grid_to_feature(before_last) |>
  model |>
  Tracker.data

grid_to_feature(before_last) |>
  model |>
  Tracker.data

function assess_model(model)
  init_grid = init_game()
  @time (init_grid, fnl, seq, cart, tf, reward) = simulate_game!(init_grid, model, 0.0)

  # obtain the discounted rewards
  y = discounted_rewards(reward, 0.99)[1:end-1];

  # score every move using  model
  updated_model = [Tracker.data(model(grid_move_to_feature(sim_seq(init_grid, seq[1:i+1], cart[1:i+1], tf[1:i+1]), seq[i])))[1] for i=1:length(seq)-1]
  display(sim_seq(init_grid, seq, cart, tf))
  g = plot(updated_model, label="model")
  plot!(y, label="discounted rewards")
  (init_grid, fnl, seq, cart, tf, reward, g)
end

ϵ = 0.5

model = nothing
@load "no_neg_model.jld2"
model

@time for i=1:10000
  global ϵ
  ϵ -= 0.0001
  simulate_one_game_and_train!(model, opt, ϵ)
  if i % 100 == 0
    @save "no_neg_model.jld2" model;
    println("epsilon: $ϵ")
    init_grid = init_game()
    res = simulate_game!(init_grid, model, 0)
    display(res[2])
    @save "no_neg_model.jld2" model;
  end
end

@time (init, fnl, seq, cart, tf, reward, g) = assess_model(model);
g

last_n_move = length(seq)-20
grid_new = sim_seq(init, seq[1:end-last_n_move], cart[1:end-last_n_move], tf[1:end - last_n_move])
grid_try = copy(grid_new)
display(grid_new)
vm = valid_moves(grid_new)
x = vec(Tracker.data(model(hcat(grid_move_to_feature.([grid_new], vm)...))))
y = x .- minimum(x)

using DataFrames
DataFrame(vm=vm, chosen=vm.== seq[end-last_n_move], x=vec(x),  prob_exp = round.(exp.(y)/sum(exp.(y)),digits=2), prob_square =((y.+1).^2)./sum((y.+1).^2))


text:
# Define functions to find ∆ω, ∆a, ∆ψ by the implicit Newton step calculation algorithm.
# Note that ω, a, ψ are fixed in this calculation.

# We can save computations by implementing the assumption ∆D = 0, and also by updating ψ
# directly without calculating ∆ψ.

export LasingSol, LasingVar

# Below (and also in nonlasing.jl), many constructors take some arguments as templates, and
# those template arguments are not directly set as fields: separate copies are created
# (usually by `similar`) and set as fields.
#
# - An alternative would be to take parameters for creating those fields and create the
# fields.  E.g., we could take the length and the element type of a vector, create a vector,
# and set it as a field.  However, we are aiming to support many different types of vectors
# (e.g., PETSc vector).  We cannot predict all the parameters needed for creating these
# vectors.  By taking an actual instance of a vector in this case, we can support any types
# of vectors.
#
# - Then, the next possibility would be to let the users to provide a copy, instead of
# automatically generating a copy inside code.  I avoid taking this route for consistency in
# regard to automatic type conversion inside constructors.  Suppose I create a custom type
# named `MyType` containing `v::Vector{Float64}` as a field.  If I instantiate `MyType` by
# giving a `vi::Vector{Int64}` as an argument, `vi` is automatically type-converted into a
# Vector{Float64} and set to `v`.  Here, the resulting `v::Vector{Float64}` is a separate
# copy from `vi::Vector{Int64}`.  On the other hand, if I instantiate `MyType` by giving a
# `vf::Vector{Float64}` as an argument, a copy is not created and the `vf` itself is set to
# `v`.  This inconsistency is confusing, because in one case changing the contents of the
# outside vector `vi` does not affect the `MyType` instance, whereas in the other case it
# does.  In order to avoid this inconsistency, I decide to always create a copy of a
# template argument.  This way, no matter whether automatic type conversion occurs or not,
# the stored field is a separate copy from the argument given to a constructor.

# Solutions to the lasing equation.
mutable struct LasingSol{VF<:AbsVecFloat,VC<:AbsVecComplexF}  # VF and VC can be PETSc vectors with N entries
    ω::VecFloat  # M real numbers: frequencies of modes (M = # of total modes, including both lasing and nonlasing)
    a²::VecFloat  # M real numbers: squared "amplitudes" of modes inside hole-burning term; this includes γ(ω)² factor
    ψ::Vector{VC}  # M complex vectors: normalized modes
    abs2ψ::Vector{VF}  # M complex vectors: absolute squares of normalized modes; this is not abs2.(ψ) because of interpolation
    iₐ::VecInt  # M integers: row indices where amplitudes are measured
    vtemp::VC  # temporary storage for N complex numbers
    wtemp::VC  # temporary storage for N complex numbers
    activated::VecBool  # activated[m] == true if mode m is just activated; used in simulate!
    active::VecBool  # active[m] == true if mode m is active (i.e., lasing)
    m_active::VecInt  # vector of active (i.e., lasing) mode indices; collection of m such that active[m] == true
    function LasingSol{VF,VC}(ω::AbsVecReal,
                              a²::AbsVecReal,
                              ψ::AbsVec{<:AbsVecNumber},
                              abs2ψ::AbsVec{<:AbsVecReal},
                              iₐ::AbsVecInteger,
                              vtemp::AbsVecComplexF) where {VF<:AbsVecFloat,VC<:AbsVecComplexF}
        M = length(ω)
        length(a²)==length(ψ)==length(abs2ψ)==length(iₐ)==M ||
            @error "length.((ω,a²,ψ,abs2ψ)) = $(length.((ω,a²,ψ,abs2ψ))) should be all same."

        N = length(vtemp)
        for m = 1:M
            length(ψ[m])==length(abs2ψ[m])==N ||
                @error "length.((ψ[$m],abs2ψ[$m],vtemp)) = $(length.((ψ[m],abs2ψ[m],vtemp))) should be all same."
        end

        return new(ω, a², ψ, abs2ψ, iₐ, vtemp, copy(vtemp), fill(false,M), fill(false,M), VecInt(undef,0))  # active[m]=false for all m: no mode is lasing
    end
end
LasingSol(ω::AbsVecReal, a²::AbsVecReal, ψ::AbsVec{VC}, abs2ψ::AbsVec{VF}, iₐ::AbsVecInteger, vtemp::VC) where {VF<:AbsVecFloat,VC<:AbsVecComplexF} =
    LasingSol{VF,VC}(ω, a², ψ, abs2ψ, iₐ, vtemp)

# To do: check if the following works for vtemp of PETSc vector type.
LasingSol(vtemp::AbsVec, M::Integer) =  # vtemp has N entries
    LasingSol(zeros(M), zeros(M), [similar(vtemp,ComplexF).=0 for m = 1:M], [similar(vtemp,Float).=0 for m = 1:M], zeros(Int,M), similar(vtemp,ComplexF))
LasingSol(N::Integer, M::Integer) = LasingSol(VecFloat(undef,N), M)

Base.length(lsol::LasingSol) = length(lsol.ω)  # number of total modes, including both lasing and nonlasing
num_active(lsol::LasingSol) = sum(lsol.active)  # number of modes currently lasing

# Note that this function changes iₐ.  Therefore, this should not be called inside the
# iteration for finding the solution for a given pump strength, because iₐ should be kept the
# same for a given pump strength in order to solve the equations with the same normalization
# conditions.
function LinearAlgebra.normalize!(lsol::LasingSol)
    for m = lsol.m_active
        ψ = lsol.ψ[m]
        iₐ = argmax(abs2, ψ)
        lsol.iₐ[m] = iₐ

        lsol.a²[m] *= abs2(ψ[iₐ])
        ψ ./= ψ[iₐ]  # make ψ[iₐ] = 1
    end
end

# Lasing equation variables that are reduced from all lasing modes
mutable struct LasingReducedVar{VF<:AbsVecFloat}  # VF can be PETSc vector with N entries
    D::Vector{VF}  # length-K vector: each entry is N-vector storing population inversion, where K is number of atomic classes
    D′::Vector{VF}  # length-K vector: each entry is N-vector storing derivative of population inversion with respect to hole-burning term
    function LasingReducedVar{VF}(D::AbsVec{<:AbsVecReal},
                                  D′::AbsVec{<:AbsVecReal}) where {VF<:AbsVecFloat}
        K = length(D)  # number of atomic classes
        length(D′)==K || @error "length(D) = $K and length(D′) = $(length(D′)) should be same"

        if K ≥ 1
            N = length(D[1])
            for k = 1:K
                length(D[k])==length(D′[k])==N ||
                    @error "length(D[$k]) = $(length(D[k])) and length(D′[$k]) = $(length(D′[k])) should be $N."
            end
        end

        return new(D, D′)
    end
end
LasingReducedVar(D::AbsVec{<:AbsVecReal}, D′::AbsVec{<:AbsVecReal}) = LasingReducedVar{eltype(D)}(D, D′)

# To do: check if the following works for vtemp of PETSc vector type.
LasingReducedVar(vtemp::AbsVec, K::Integer) = LasingReducedVar([similar(vtemp,Float).=0 for k=1:K], [similar(vtemp,Float).=0 for k=1:K])  # vtemp has N entries
LasingReducedVar(N::Integer, K::Integer) =  LasingReducedVar(VecFloat(undef,N), K)

# Initialize the variables that are NOT specific to a specific lasing mode.  These variables
# are constructed by summing up the contributions from all lasing modes.
function init_reduced_var!(rvar::LasingReducedVar, lsol::LasingSol, gp::GainProfile)
    # abs2ψ should be updated before calling this function; see init_lvar!.

    hb = lsol.vtemp  # temporary storage for hole-burning term
    for k = 1:length(gp)  # atomic classes
        # Calculate the (1 + hole-burning term) for atomic class k.
        abs2γ = gp.abs2gain[k]  # |γ²| function for atomic class k
        hole_burning!(hb, abs2γ, lsol.ω, lsol.a², lsol.abs2ψ)

        # Calculate the population inversion and the derivative of the population inversion
        # for the class-k atoms.
        #
        # hb is a complex vector storing real, and it is faster to take its real part first
        # and perform the further calculation on it.
        wt = gp.wt[k]  # weight function for atomic class k
        rvar.D[k] .= wt.(gp.D₀) ./ real.(hb)
        rvar.D′[k] .= (-1.0) .* wt.(gp.D₀) ./ abs2.(real.(hb))
    end

    return nothing
end

# Lasing equation variables that are unique to each lasing mode
mutable struct LasingModalVar{MXS<:AbstractMaxwellSolver,VC<:AbsVecComplexF}  # VC can be PETSc vector
    # Consider adding εeff as a field, in case the user wants to solve the "linear" eigenvalue
    # equation for some reason.
    mxsolver::MXS
    ∂f∂ω::VC
end

# To do: check if the following works for vtemp of PETSc vector type.
LasingModalVar(mxsolver_temp::AbstractMaxwellSolver, vtemp::AbsVec) =  # for, e.g., PETSc vector vtemp; vtemp has N entries
    LasingModalVar(similar(mxsolver_temp), similar(vtemp,ComplexF))

LasingModalVar(mxsolver_temp::AbstractMaxwellSolver, N::Integer) =
    LasingModalVar(similar(mxsolver_temp), VecComplexF(undef,N))

# Initialize the variables that are specific to a specific lasing mode.
function init_modal_var!(mvar::LasingModalVar,
                         m::Integer,  # index of lasing mode
                         lsol::LasingSol,
                         D::AbsVec{<:AbsVecFloat},
                         gp::GainProfile,
                         εc::AbsVecComplexF)
    init_modal_var_impl!(mvar.mxsolver, mvar.∂f∂ω, lsol.vtemp, lsol.ω[m], lsol.ψ[m], D, gp, εc)
end

# This is used in nonlasing.jl as well.
function init_modal_var_impl!(mxsolver::AbstractMaxwellSolver,
                              ∂f∂ω::AbsVecComplexF,
                              vtemp::AbsVecComplexF,
                              ω::Number,
                              ψ::AbsVecComplexF,
                              D::AbsVec{<:AbsVecFloat},
                              gp::GainProfile,
                              εc::AbsVecComplexF)
    K = length(gp)  # number of atomic classes
    γ = gp.gain  # K-vector whose kth entry is γ function of atomic class k
    γ′ = gp.gain′  # K-vector whose kth entry is γ′ function of atomic class k

    # Calculate the effective permittivity.
    εeff = vtemp  # temporary storage for effective permitivity: εc + ∑ₖ(γₖ Dₖ)
    εeff .= εc
    for k = 1:K
        εeff .+= γ[k](ω) .* D[k]
    end

    init_mxsolver!(mxsolver, ω, εeff)

    # Calculate ∂f∂ω: derivative of SALT function w.r.t. ω, ignoring ω-dependence of D.
    ∂f∂ω .= 2ω .* εeff
    for k = 1:K
        σ = ω^2 * γ′[k](ω)  # scalar
        ∂f∂ω .+= σ .* D[k]
    end
    ∂f∂ω .*= ψ

    return nothing
end

mutable struct LasingConstraint
    A::MatFloat  # constraint matrix
    b::VecFloat  # constraint vector
    ∆ωa²::VecFloat  # vector of calculated ∆ω and ∆a²: ∆ω[m] = ∆ωa²[2m-1] and ∆a²[m] = ∆ωa²[2m]
    LasingConstraint(M::Integer) = new(MatFloat(undef,2M,2M), VecFloat(undef,2M), VecFloat(undef,2M))
end

# Components necessary for fixed-point iteration
mutable struct LasingVar{MXS<:AbstractMaxwellSolver,VF<:AbsVecFloat,VC<:AbsVecComplexF}
    rvar::LasingReducedVar{VF}
    mvar_vec::Vector{LasingModalVar{MXS,VC}}
    cst::LasingConstraint
    inited::Bool
end
LasingVar(mxsolver_temp::AbstractMaxwellSolver, vtemp::AbsVec, K::Integer, M::Integer) =  # K: number of atomic classes, M: number of modes
    LasingVar(LasingReducedVar(vtemp, K), [LasingModalVar(mxsolver_temp, vtemp) for m = 1:M], LasingConstraint(M), false)
LasingVar(mxsolver_temp::AbstractMaxwellSolver, N::Integer, K::Integer, M::Integer) = LasingVar(mxsolver_temp, VecFloat(undef,N), K, M)

function init_lvar!(lvar::LasingVar, lsol::LasingSol, gp::GainProfile, εc::AbsVecComplexF)
    rvar = lvar.rvar
    mvar_vec = lvar.mvar_vec

    # Update abs2ψ.
    for m = lsol.m_active
        abs2_interp!(lsol.abs2ψ[m], mvar_vec[m].mxsolver, lsol.ψ[m])
        # lsol.abs2ψ[m] .= abs2.(lsol.ψ[m])
    end

    init_reduced_var!(rvar, lsol, gp)

    for m = lsol.m_active
        init_modal_var!(mvar_vec[m], m, lsol, rvar.D, gp, εc)
    end

    lvar.inited = true

    return nothing
end

function norm_leq_impl(lsol::LasingSol, mvar_vec::AbsVec{<:LasingModalVar})
    leq = 0.0
    b = lsol.vtemp
    for m = lsol.m_active
        mxsolver = mvar_vec[m].mxsolver
        ψ = lsol.ψ[m]
        mxapply!(b, mxsolver, ψ)  # b = A * ψ
        leq = max(leq, norm(b))  # 2-norm for each mode, ∞-norm between modes
    end

    return leq  # return 0.0 if lsol.m_active is empty
end

function norm_leq(lsol::LasingSol, lvar::LasingVar, gp::GainProfile)
    lvar.inited || @error "lvar is uninitialized: call init_lvar!(...) first."

    return norm_leq_impl(lsol, lvar.mvar_vec)
end

# Create the mth constraint equation for ∆ω and ∆a.
function set_constraint!(cst::LasingConstraint,
                         lsol::LasingSol,
                         indₘ::Integer,  # index of m in m_active
                         mvar::LasingModalVar,  # modal variables for mth lasing mode
                         rvar::LasingReducedVar,
                         gp::GainProfile)
    # Retrieve necessary variables for constructing the constraint.
    m = lsol.m_active[indₘ]  # modal index of equation of interest
    L = num_active(lsol)  # number of currently lasing modes
    K = length(gp)  # number of atomic classes

    A = cst.A  # left-hand-side matrix of constraint equation
    b = cst.b  # right-hand-side vector of constraint equation

    ω = lsol.ω  # M-vector whose mth entry is current approximate mth eigenfrequency
    a² = lsol.a²  # M-vector whose mth entry is squared amplitude of current approximate mth eigenmode
    ψ = lsol.ψ  # M-vector whose mth entry is current approximate mth eigenmode
    abs2ψ = lsol.abs2ψ  # M-vector whose mth entry is absolute square of current approximate mth eigenmode

    ωₘ = ω[m]
    ωₘ² = ωₘ^2
    ψₘ = ψ[m]
    iₐₘ = lsol.iₐ[m]  # where mth eigenmode is normalized

    vtemp = lsol.vtemp  # scratch space 1
    wtemp = lsol.wtemp  # scratch space 2

    γ = gp.gain  # K-vector whose kth entry is γ function of atomic class k
    abs2γ = gp.abs2gain  # K-vector whose kth entry is |γ|² function of atomic class k
    abs2γ′ = gp.abs2gain′  # K-vector whose kth entry is ∂|γ|²/∂ω function of atomic class k

    D′ = rvar.D′  # K-vector whose kth entry is derivative of population inversion of class-k atoms w.r.t. hole-burning term

    # Calculate the iₐth row of A⁻¹ and keep it as a column vector
    eᵢₐ = vtemp
    eᵢₐ .= 0
    eᵢₐ[iₐₘ] = 1
    r = wtemp  # storage for row vector
    mxsolve_transpose!(r, mvar.mxsolver, eᵢₐ)

    # Construct A and b.  Note that they are initialized to zero outside the present
    # function, by update_lsol_impl!.

    # Set the mth complex row of right-hand-side vector b of the constraint.
    @assert ψₘ[iₐₘ] ≈ 1
    b[2indₘ-1] = 1.0  # real(ψₘ[iₐₘ])
    b[2indₘ] = 0.0  # imag(ψₘ[iₐₘ])

    # Set the mth complex row of the left-hand-side matrix A of the constraint.
    for indⱼ = 1:L
        j = lsol.m_active[indⱼ]  # jth complex column

        ωⱼ = ω[j]
        a²ⱼ = a²[j]
        abs2ψⱼ = abs2ψ[j]

        # Set the coefficient of ∆ωⱼ.
        vtemp .= 0
        ωₘ²a²ⱼ = ωₘ² * a²ⱼ  # scalar
        for k = 1:K
            σ = (ωₘ²a²ⱼ * abs2γ′[k](ωⱼ)) * γ[k](ωₘ)  # scalar: (real * real) * complex
            vtemp .+= σ .* D′[k]
        end
        vtemp .*= ψₘ .* abs2ψⱼ

        # The mth complex row has an additional contribution in the mth complex column from
        # the ωₘ-dependence outside the ω-dependent D.
        j==m && (vtemp .+= mvar.∂f∂ω)

        mω = BLAS.dotu(r, vtemp)
        A[2indₘ-1, 2indⱼ-1] = real(mω)
        A[2indₘ, 2indⱼ-1] = imag(mω)

        # Set the coefficient of ∆a²ⱼ.
        vtemp .= 0
        for k = 1:K
            σ = (ωₘ² * abs2γ[k](ωⱼ)) * γ[k](ωₘ)  # scalar: (real * real) * complex
            vtemp .+= σ .* D′[k]
        end
        vtemp .*= ψₘ .* abs2ψⱼ

        ma² = BLAS.dotu(r, vtemp)
        A[2indₘ-1, 2indⱼ] = real(ma²)
        A[2indₘ, 2indⱼ] = imag(ma²)
    end

    return nothing
end

# Update ψ for the mth mode.  This function does not update ω and a².
function update_ψₘ!(lsol::LasingSol,
                    ∆ωa²::VecFloat,
                    indₘ::Integer,  # index of m in m_active
                    mvar::LasingModalVar,
                    rvar::LasingReducedVar,
                    gp::GainProfile)
    # Retrieve necessary variables for constructing the constraint.
    m = lsol.m_active[indₘ]  # modal index of equation of interest
    L = num_active(lsol)  # number of currently lasing modes
    K = length(gp)  # number of atomic classes

    ω = lsol.ω  # M-vector whose mth entry is current approximate mth eigenfrequency
    a² = lsol.a²  # M-vector whose mth entry is squared amplitude of current approximate mth eigenmode
    ψ = lsol.ψ  # M-vector whose mth entry is current approximate mth eigenmode
    abs2ψ = lsol.abs2ψ  # M-vector whose mth entry is absolute square of current approximate mth eigenmode

    ωₘ = ω[m]
    ωₘ² = ωₘ^2
    ψₘ = ψ[m]
    iₐₘ = lsol.iₐ[m]  # where mth eigenmode is normalized

    vtemp = lsol.vtemp  # scratch space 1
    wtemp = lsol.wtemp  # scratch space 1

    γ = gp.gain  # K-vector whose kth entry is γ function of atomic class k
    abs2γ = gp.abs2gain  # K-vector whose kth entry is |γ|² function of atomic class k
    abs2γ′ = gp.abs2gain′  # K-vector whose kth entry is ∂|γ|²/∂ω function of atomic class k

    D′ = rvar.D′  # K-vector whose kth entry is derivative of population inversion of class-k atoms w.r.t. hole-burning term

    wtemp .= 0  # initialize for mth row
    for indⱼ = 1:L
        j = lsol.m_active[indⱼ]  # jth complex column

        ωⱼ = ω[j]
        a²ⱼ = a²[j]
        abs2ψⱼ = abs2ψ[j]

        ∆ωⱼ = ∆ωa²[2indⱼ-1]
        ∆a²ⱼ = ∆ωa²[2indⱼ]

        # Add the contribution of ∆ωⱼ.
        vtemp .= 0  # initialize for ∆ωⱼ
        ωₘ²a²ⱼ = ωₘ² * a²ⱼ  # scalar
        for k = 1:K
            σ = (∆ωⱼ * ωₘ²a²ⱼ * abs2γ′[k](ωⱼ)) * γ[k](ωₘ)  # scalar: (real * real * real) * complex
            vtemp .+= σ .* D′[k]
        end
        vtemp .*= abs2ψⱼ  # ψₘ will be multiplied outside for loop, because it is common for both ∆ω and ∆a²
        wtemp .+= vtemp

        # Add the contribution of ∆a²ⱼ.
        vtemp .= 0  # initialize for ∆a²ⱼ
        for k = 1:K
            σ = (∆a²ⱼ * ωₘ² * abs2γ[k](ωⱼ)) * γ[k](ωₘ)  # scalar: (real * real * real) * complex
            vtemp .+= σ .* D′[k]
        end
        vtemp .*= abs2ψⱼ  # ψₘ will be multiplied outside for loop, because it is common for both ∆ω and ∆a²
        wtemp .+= vtemp
    end
    wtemp .*= ψₘ

    # The mth complex row has an additional contribution from the ωₘ-dependence outside the
    # ω-dependent D.
    ∆ωₘ = ∆ωa²[2indₘ-1]
    wtemp .+= ∆ωₘ .* mvar.∂f∂ω

    # @info "‖A‖₁ = $(opnorm(mvar.A,1)), ‖vtemp‖ = $(norm(vtemp))"
    mxsolve!(ψₘ, mvar.mxsolver, wtemp)  # ψₘ .= mvar.mxsolver \ wtemp

    # Normalize ψₘ just in case ψₘ[iₐₘ] ≠ 1.
    # ψₘ[iₐₘ]≈1 || @warn "lasing mode m = $m is slightly nonnormal: |ψₘ[iₐₘ]-1| = $(abs(ψₘ[iₐₘ]-1)).  The mode will be renormalized."
    # @info "ψₘ[iₐₘ] = $(ψₘ[iₐₘ])"
    ψₘ ./= ψₘ[iₐₘ]

    return nothing
end

# lvar should be already initialized by init_lvar! before starting the fixed-point iteration.
function update_lsol!(lsol::LasingSol, lvar::LasingVar, gp::GainProfile)
    lvar.inited || @error "lvar is uninitialized: call init_lvar!(...) first."
    update_lsol_impl!(lsol, lvar.mvar_vec, lvar.rvar, lvar.cst, gp)
    lvar.inited = false

    return nothing
end

# Fixed-point equation for lsol.  Calculate ∆ω and ∆a, and then ∆ψ, and move lsol by them.
#
# This function is expensive to call, because it involves 2Mₗ linear solves, where Mₗ is the
# number of lasing modes.  Therefore, try not to call this function unnecessarily.
#
# init_lvar! should be called before using this function to make lvar prepared.  However,
# init_lvar! is not exported in order to force checking the norm by norm_leq to avoid
# calling this function when the norm is small enough.
function update_lsol_impl!(lsol::LasingSol,
                           mvar_vec::AbsVec{<:LasingModalVar},  # should be already initialized
                           rvar::LasingReducedVar,  # should be already initialized
                           cst::LasingConstraint,
                           gp::GainProfile)
    # Construct the constraint equation on ∆ω and ∆a.
    L = num_active(lsol)

    A = cst.A
    b = cst.b
    ∆ωa² = cst.∆ωa²

    A .= 0
    b .= 0
    ∆ωa² .= 0

    # Below, we iterate over the contiguous indices indₘ to retrieve the currently lasing
    # mode index m from m_active, instead of directly iterating over m_active.  This is
    # because we cannot use a submatrix @view(A[ind,ind]) in ldiv! when ind is not a
    # contiguous range of indices.  See https://github.com/JuliaLang/julia/issues/30097.
    for indₘ = 1:L
        m = lsol.m_active[indₘ]
        set_constraint!(cst, lsol, indₘ, mvar_vec[m], rvar, gp)
    end

    # Calculate ∆ω and ∆a.
    # @info "A = $(A[1:2L,1:2L]), b = $(b[1:2L]), ∆ωa² = $(∆ωa²[1:2L])"
    LU = lu!(@view(A[1:2L,1:2L]))
    ldiv!(@view(∆ωa²[1:2L]), LU, @view(b[1:2L]))

    # Update ψ.
    # @info "lsol.ω = $(lsol.ω), lsol.a² = $(lsol.a²), lsol.m_active = $(lsol.m_active)"
    for indₘ = 1:L
        m = lsol.m_active[indₘ]
        # @info "before apply: ‖lsol.ψ[$m]‖ = $(norm(lsol.ψ[m]))"
        update_ψₘ!(lsol, ∆ωa², indₘ, mvar_vec[m], rvar, gp)
        # @info "after apply: ‖lsol.ψ[$m]‖ = $(norm(lsol.ψ[m]))"
    end

    # Update ω and a².
    #
    # abs2ψ is not updated here.  In each iteration, the Anderson acceleration evaluates the
    # fixed-point function g once to get g(xₖ), but this g(xₖ) is not the new solution xₖ₊₁.
    # xₖ₊₁ is formed as an affine combinination of g(xₖ), ..., g(xₖ₋ₘ).  Therefore,
    # abs2ψ evaluated at the end of the evaluation of the fixed point function g is abs2ψ of
    # g(xₖ), not the abs2ψ of xₖ₊₁.  Therefore, in order to calculate abs2ψ for xₖ₊₁, which
    # is used inside g for evaluating g(xₖ₊₁), abs2ψ should be updated at the beginning of the
    # fixed-point function (where abs2ψ is used for the first time in that iteration step)
    # rather than at the end.  This is inside init_lvar!, whose evaluation in checked at the
    # beginning of update_lsol!.
    #
    # The following loop should not be merged with the loop above, because update_ψₘ! uses
    # un-updated ω and a².
    for indₘ = 1:L
        m = lsol.m_active[indₘ]

        ∆ωₘ = ∆ωa²[2indₘ-1]
        ∆a²ₘ = ∆ωa²[2indₘ]

        lsol.ω[m] += ∆ωₘ
        lsol.a²[m] += ∆a²ₘ
    end

    return nothing
end


text:
"""
Materialize a `Projection` into a `DataFrame` by generating each output column
in sequence.
"""
function materialize(node::Projection)
    out = DataFrames.DataFrame()
    # TODO: Perfomanance optimization
    #    If there are no cross-dependencies, generate the output in parallel.
    for e in node.expressions
        out[e.alias] = generate_column(node.source, e)
    end
    out
end