text:
function initialize_dict(r,c,n)
    # initialize memory addr and memory neighbor addr
    mem_dict = mem_layout(r,c,n)
    dict = Dict{String,Any}()
    # move to a bigger dict
    dict["memory_dict"] = mem_dict
    dict["input_handler"] = 1
    return dict
end

function mem_layout(row,col,count)
    row_spacing = floor((row-1)/(count/4)) #leave one row for input handler
    col_spacing = floor(col/(count/4))
    row_addr = Int64[]
    col_addr = Int64[]
    for i = 0:(count/4)-1
        push!(row_addr,(i*row_spacing)+4)
    end
    for i = 0:(count/4)-1
        push!(col_addr,(i*col_spacing)+4)
    end
    mem_dict = Dict{Mapper2.CartesianIndex{3},Array{Mapper2.CartesianIndex{3},1}}()
    for r in row_addr
        mem_addr = CartesianIndex(r,1,1)
        memproc_addr = [CartesianIndex(r,2,1),CartesianIndex(r+1,2,1)]
        mem_dict[mem_addr] = memproc_addr
        mem_addr = CartesianIndex(r,col+2,1)
        memproc_addr = [CartesianIndex(r,col+1,1),CartesianIndex(r+1,col+1,1)]
        mem_dict[mem_addr] = memproc_addr
    end
    for c in col_addr
        mem_addr = CartesianIndex(1,c,1)
        memproc_addr = [CartesianIndex(2,c,1),CartesianIndex(2,c+1,1)]
        mem_dict[mem_addr] = memproc_addr
        mem_addr = CartesianIndex(row+2,c,1)
        memproc_addr = [CartesianIndex(row+1,c,1),CartesianIndex(row+1,c+1,1)]
        mem_dict[mem_addr] = memproc_addr
    end

    return mem_dict
end

function generic(row::Int64,
                 col::Int64,
                 lev::Int64,
                 num_memories::Int64,
                 A = KCStandard,
                 num_links = 2)

    dict = initialize_dict(row,col,num_memories)

    # check for memory_dict key in the bigger dictionary
    if haskey(dict, "input_handler")
        num_in = dict["input_handler"]
    end

    if haskey(dict, "memory_dict")
        mem_dict = dict["memory_dict"]
    else
        # if dict does not contain memory_dict, initialize an empty array
        mem_dict = Dict{Mapper2.CartesianIndex{3},Vector{Mapper2.CartesianIndex{3}}}()
    end

    # Check the dimension
    if lev == 1
        dim = 2
    elseif lev > 1
        dim = 3
    else
        error("The level argument is invalid.")
    end

    # direction definitions
    directions2d = ("north", "east", "south", "west")
    directions3d = ("north", "east", "south", "west", "top", "bottom")
    directions = (directions2d,directions3d)

    # Start with a new component
    arch = TopLevel{A,3}("generic")

    # make empty arrays for memory addresses and memory neighbor addresses
    mem_array = CartesianIndex[]
    mem_proc_array = CartesianIndex[]
    # move the addresses from mem_dict to arrays
    for (key,value) in mem_dict
        push!(mem_array,key)
        for addr in value
            push!(mem_proc_array,addr)
        end
    end

	####################
    # Normal Processor #
    ####################
    # Get a processor tile to instantiate
	processor = build_processor_tile_generic(dim,directions)
    # Instantiate it at the required addresses
  	for r in 2:row+1, c in 2:col+1, l in 1:lev
        if CartesianIndex(r,c,l) in mem_proc_array || CartesianIndex(r,c,l) in mem_array
            continue # avoid the addresses in mem_dict
        end
		add_child(arch, processor, CartesianIndex(r,c,l))
  	end

	####################
   	# Memory Processor #
   	####################
   	memory_processor = build_memory_processor_tile_generic(dim,directions)
   	for addr in mem_proc_array
	  	add_child(arch, memory_processor, addr)
   	end

    #################
    # 2 Port Memory #
    #################
    memory_2port = build_memory_2port()
    for addr in mem_array
        add_child(arch, memory_2port, addr)
    end

    #################
    # Input Handler #
    #################
    input_handler = build_input_handler(num_links)
    s = floor(row/num_in) # spacing
    for i = 0:num_in-1
        row = Int(2+(s*i))
        println(row)
        add_child(arch, input_handler, CartesianIndex(row,1,1))
    end
    ##################
    # Output Handler #
    ##################
    output_handler = build_output_handler(num_links)
    add_child(arch, output_handler, CartesianIndex(2,col+2,1))

    #######################
    # Global Interconnect #
    #######################
    connect_processors_generic(arch,dim)
    connect_memories_generic(arch)
    return arch

end

##############################
#        PROCESSOR
##############################
function build_processor_generic(dimension::Int64,
directions::Tuple{NTuple{4,String},NTuple{6,String}})

    # Build the metadata dictionary for the processor component
    metadata = Dict{String,Any}()
    metadata["attributes"] = ["processor"]
    component = Component(  "standard_processor",
                            primitive = "", metadata = metadata)
    # Add the input fifos
    add_port(component, "fifo", "input", 2)
    # Add the output ports
    for dir in directions[dimension-1]
        add_port(component, dir, "output", 2)
    end
    # Add the dynamic circuit switched network
    add_port(component, "dynamic", "output")
    # Add memory ports. Will only be connected in the memory processor tile.
    add_port(component, "memory_in", "input")
    add_port(component, "memory_out", "output")
    # Return the created type
    return component

end

##################
# COMPLEX BLOCKS #
##################
function build_processor_tile_generic(dimension::Int64,
directions::Tuple{NTuple{4,String},NTuple{6,String}})
    # Working towards parameterizing this. For now, just leave this at two
    # because the "processor" components aren't parameterized for the number
    # of ports. This should be easy to fix though.
    num_links = 2
    # Create a new component for the processor tile
    # No need to set the primtiive class or metadata because we won't
    # be needing it.
    comp = Component("processor_tile")
    # Add the circuit switched ports
    for dir in directions[dimension-1]
        for (suffix,class)  in zip(("_in", "_out"), ("input", "output"))
            port_name = join((dir, suffix))
            add_port(comp, port_name, class, num_links)
        end
    end
    # Add memory ports - only memory processor tiles will have the necessary
    # "memory_processor" attribute in the core to allow memory applicationa
    # to be mapped to them.
    add_port(comp, "memory_in", "input")
    add_port(comp, "memory_out", "output")
    # Instantiate the processor primitive
    add_child(comp, build_processor_generic(dimension,directions), "processor")
    # Instantiate the directional routing muxes
    routing_mux = build_mux(length(directions[dimension-1]),1)
    for dir in directions[dimension-1]
        name = join((dir, "_mux"))
        add_child(comp, routing_mux, name, num_links)
    end
    # Instantiate the muxes routing data to the fifos
    add_child(comp, build_mux(  (length(directions[dimension-1])*2)+1,1),
                                "fifo_mux", 2)

    # Interconnect - Don't attach metadata and let the routing routine fill in
    # defaults to intra-tile routing.
    add_link(comp, "processor.memory_out", "memory_out")
    add_link(comp, "memory_in", "processor.memory_in")

    # Connect outputs of muxes to the tile outputs
    for dir in directions[dimension-1], i = 0:num_links-1
        mux_port = join((dir, "_mux[",string(i),"].out[0]"))
        tile_port = join((dir, "_out[",string(i),"]"))
        add_link(comp, mux_port, tile_port)
    end
    # Circuit switch output links
    for dir in directions[dimension-1], i = 0:num_links-1
        processor_port = join(("processor.",dir,"[",string(i),"]"))
        mux_port = join((dir, "_mux[",string(i),"].in[0]"))
        add_link(comp,processor_port,mux_port)
    end

    # Input Links
    add_link(comp, "fifo_mux[0].out[0]", ["processor.fifo[0]"])
    add_link(comp, "fifo_mux[1].out[0]", ["processor.fifo[1]"])
    # Connect input ports to inputs of muxes
    dir_count = 0
    fifo_count = 0
    dir_tracker = Dict(d => 0 for d in directions[dimension-1])
    for dir in directions[dimension-1]
        dir_count += 1
        for i = 0:num_links-1
            fifo_count += 1
            tile_port = join((dir, "_in[",string(i),"]"))
            sink_ports = String[]
            for d in directions[dimension-1]
                if d == dir
                    dir_tracker[d] = 1
                    continue
                end
                mux_port = join((d, "_mux[",string(i),"].in[",
                                string(dir_count-dir_tracker[d]),"]"))
                push!(sink_ports,mux_port)
            end
            push!(sink_ports,join(("fifo_mux[0].in[",string(fifo_count-1),"]")))
            push!(sink_ports,join(("fifo_mux[1].in[",string(fifo_count-1),"]")))
            add_link(comp, tile_port, sink_ports)
        end
    end
    return comp

end

function build_memory_processor_tile_generic(dimension::Int64,
directions::Tuple{NTuple{4,String},NTuple{6,String}})
    # Get a normal processor and add the memory ports to it.
    tile = build_processor_tile_generic(dimension,directions)
    # Need to add the memory processor attribute the processor.
    push!(tile.children["processor"].metadata["attributes"], "memory_processor")
    return tile
end

function connect_processors_generic(tl,dimension)
    # General rule - we're looking for the attribute "processor" to be somewhere
    # in the component stack. If so, we'll try to connect all of the circuit
    # switched ports.
    vals = ["processor", "input_handler", "output_handler"]
    fn = x -> search_metadata!(x, "attributes", vals, oneofin)
    src_rule = fn
    dst_rule = fn

    # Build metadata dictionary for capacity and cost
    metadata = Dict(
        "cost"      => 1.0,
        "capacity"  => 1,
        "network"   => ["circuit_switched"]
    )

    # Create offset rules.
    # Offsets are just unit steps in four directions.
    offsets2d = [CartesianIndex(-1,0,0),
                 CartesianIndex(1,0,0),
                 CartesianIndex(0,1,0),
                 CartesianIndex(0,-1,0)
                ]

    offsets3d = [CartesianIndex(-1,0,0),
                 CartesianIndex(1,0,0),
                 CartesianIndex(0,1,0),
                 CartesianIndex(0,-1,0),
                 CartesianIndex(0,0,1),
                 CartesianIndex(0,0,-1)
                ]

    offsets = (offsets2d,offsets3d)

    #=
    Create two tuples for the source ports and destination ports. In general,
    if the source link is going out of the north port, the destionation will
    be coming in the south port.
    =#

    src_dirs2d = ("north", "south", "east", "west")
    dst_dirs2d = ("south", "north", "west", "east")
    src_dirs3d = ("north", "south", "east", "west", "top", "bottom")
    dst_dirs3d = ("south", "north", "west", "east", "bottom", "top")
    src_dirs = (src_dirs2d,src_dirs3d)
    dst_dirs = (dst_dirs2d,dst_dirs3d)

    offset_rules = []
    for (offset, src, dst) in zip(offsets[dimension-1], src_dirs[dimension-1],
    dst_dirs[dimension-1])
        src_ports = String[]
        dst_ports = String[]
        # Iterate through the number of source and destination ports.
        for i in 0:1
            src_port = "$(src)_out[$i]"
            push!(src_ports, src_port)
            dst_port = "$(dst)_in[$i]"
            push!(dst_ports, dst_port)
        end
        # Create the offset rule and add it to the collection
        orule = Splatter(Iterators.product( (offset,),
                                            zip(src_ports,dst_ports)
                                          ))
        connection_rule(tl, orule, src_rule, dst_rule, metadata = metadata)
    end
    # Create offset rules for the input and output handlers.
    # Input and output handlers only appear on the left and right hand sides
    # of the array, so only need the "east" and "west" directions.
    src_dirs = ("east","west")
    dst_dirs = ("west","east")
    # Links can go both directions, so make the offsets an array
    offsets = [CartesianIndex(0,1,0), CartesianIndex(0,-1,0)]
    for (offset, src, dst) in zip(offsets, src_dirs, dst_dirs)
        src_ports = String[]
        dst_ports = String[]
        for i in 0:1
            # Input handler -> processor
            src_port = "out[$i]"
            dst_port = "$(dst)_in[$i]"
            push!(src_ports, src_port)
            push!(dst_ports, dst_port)
            # processor -> output handler
            src_port = "$(src)_out[$i]"
            dst_port = "in[$i]"
            push!(src_ports, src_port)
            push!(dst_ports, dst_port)
        end
        orule = Splatter(Iterators.product( (offset,),
                                            zip(src_ports,dst_ports)
                                          ))

        connection_rule(tl, orule, src_rule, dst_rule, metadata = metadata)
    end
    return nothing
end

function connect_memories_generic(tl)
    # Create metadata dictionary for the memory links.
    metadata = Dict(
        "cost"      => 1.0,
        "capacity"  => 1,
        "network"   => ["memory"],
   )
    ###########################
    # Connect 2 port memories #
    ###########################
    proc_rule = x -> search_metadata!(x, "attributes", "memory_processor", in)
    mem_rule  = x -> search_metadata!(x, "attributes", "memory_2port", in)
    # Make connections from memory to memory-processors
    offset_rules = [
        (CartesianIndex(-1, 0,0), "out[0]", "memory_in"),
        (CartesianIndex(-1, 1,0), "out[1]", "memory_in"),
        (CartesianIndex( 1, 0,0), "out[0]", "memory_in"),
        (CartesianIndex( 1, 1,0), "out[1]", "memory_in"),
        (CartesianIndex( 0, 1,0), "out[0]", "memory_in"),
        (CartesianIndex( 1, 1,0), "out[1]", "memory_in"),
        (CartesianIndex( 0,-1,0), "out[0]", "memory_in"),
        (CartesianIndex( 1,-1,0), "out[1]", "memory_in"),
       ]
    connection_rule(tl, offset_rules, mem_rule, proc_rule, metadata = metadata)

    # Make connections from memory-processors to memories.
    offset_rules = [
        (CartesianIndex( 1, 0,0), "memory_out", "in[0]"),
        (CartesianIndex( 1,-1,0), "memory_out", "in[1]"),
        (CartesianIndex(-1, 0,0), "memory_out", "in[0]"),
        (CartesianIndex(-1,-1,0), "memory_out", "in[1]"),
        (CartesianIndex( 0,-1,0), "memory_out", "in[0]"),
        (CartesianIndex(-1,-1,0), "memory_out", "in[1]"),
        (CartesianIndex( 0, 1,0), "memory_out", "in[0]"),
        (CartesianIndex(-1, 1,0), "memory_out", "in[1]"),
       ]
    connection_rule(tl, offset_rules, proc_rule, mem_rule, metadata = metadata)

    return nothing
end


text:
struct Astrometry{N,T<:Number}
    epoch::SVector{N,T}
    ra::SVector{N,T}
    dec::SVector{N,T}
    σ_ra::SVector{N,T}
    σ_dec::SVector{N,T}
end
export Astrometry

function Astrometry(observations::NamedTuple...)
    if any(<(1), abs.(getproperty.(observations, :ra))) ||
        any(<(1), abs.(getproperty.(observations, :dec)))
        @warn "Confirm that astrometry is entered in milliarcseconds (very small values detected)"
    end
    return Astrometry(
        SVector(getproperty.(observations, :epoch)),
        SVector(getproperty.(observations, :ra)),
        SVector(getproperty.(observations, :dec)),
        SVector(getproperty.(observations, :σ_ra)),
        SVector(getproperty.(observations, :σ_dec)),
    )
end

struct Photometry{N,T<:Number}
    band::SVector{N,Symbol}
    phot::SVector{N,T}
    σ_phot::SVector{N,T}
end
export Photometry

function Photometry(observations::NamedTuple...)
    return Photometry(
        SVector(getproperty.(observations, :band)),
        SVector(getproperty.(observations, :phot)),
        SVector(getproperty.(observations, :σ_phot)),
    )
end

struct ProperMotionAnom{N,T<:Number}
    ra_epoch::SVector{N,T}
    dec_epoch::SVector{N,T}
    pm_ra::SVector{N,T}
    pm_dec::SVector{N,T}
    σ_pm_ra::SVector{N,T}
    σ_pm_dec::SVector{N,T}
end
export ProperMotionAnom
function ProperMotionAnom(observations::NamedTuple...)
    T = promote_type(typeof.(values(first(observations)))...)
    return ProperMotionAnom{length(observations),T}(
        SVector(getproperty.(observations, :ra_epoch)),
        SVector(getproperty.(observations, :dec_epoch)),
        SVector(getproperty.(observations, :pm_ra)),
        SVector(getproperty.(observations, :pm_dec)),
        SVector(getproperty.(observations, :σ_pm_ra)),
        SVector(getproperty.(observations, :σ_pm_dec)),
    )
end

"""
    Images(...)

A block of images of a system. Pass a vector of named tuples with the following fields:
* band
* image
* epoch
* platescale
* contrast

For example:
```julia
Images(
    (; epoch=1234.0, band=:J, image=readfits("abc.fits"), platescale=19.4)
)
```
Contrast can be a function that returns the 1 sigma contrast of the image from a separation in mas to the same units as the image file.
Or, simply leave it out and it will be calculated for you.
Epoch is in MJD.
Band is a symbol which matches the one used in the planet's `Priors()` block.
Platescale is in mas/px.
"""
struct Images{TImg,TCont}
    band::Vector{Symbol}
    image::Vector{TImg}
    epoch::Vector{Float64}
    platescale::Vector{Float64}
    contrast::Vector{TCont}
end
export Images
function Images(observations::NamedTuple...)
    observations_prepared = map(observations) do obs
        if !hasproperty(obs, :contrast)
            @info "Measuring contrast from image"
            contrast = contrast_interp(obs.image)
            obs = merge(obs, (;contrast))
        end
        return obs
    end
    return Images(
        collect(getproperty.(observations_prepared, :band)),
        collect(getproperty.(observations_prepared, :image)),
        collect(getproperty.(observations_prepared, :epoch)),
        collect(getproperty.(observations_prepared, :platescale)),
        collect(getproperty.(observations_prepared, :contrast)),
    )
end

"""
    Priors(key=dist, ...)

A group of zero or more priors passed by keyword arguments.
The priors must be univariate distributions from the Distributions.jl
package.
"""
struct Priors{T}
    priors::T
end
export Priors
# Basically just wrap a named tuple
Priors(;priors...) = Priors{typeof(priors)}(priors)

function Base.show(io::IO, mime::MIME"text/plain", @nospecialize priors::Priors)
    println(io, "Priors:")
    for (k,prior) in zip(keys(priors.priors), values(priors.priors))
        print(io, "  ", k, "\t")
        show(io, mime, prior)
        print(io, "\n")
    end
end

"""
    Derived(key=func, ...)

A group of zero or more functions that resolve to a parameter for the model.
Derived parameters must be functions accepting one argument if applied to a System,
or two arguments if applied to a Planet.
Must always return the same value for the same input (be a pure function) and
support autodifferentiation.
"""
struct Derived{T}
    variables::T
end
const Deterministic = Derived
export Derived, Deterministic
function Derived(;variables...)
    # Basically just wrap a named tuple
    return Derived(
        NamedTuple(variables),
    )
end
function Base.show(io::IO, mime::MIME"text/plain", @nospecialize det::Derived)
    print(io, "Derived:\n  ")
    for k in keys(det.variables)
        print(io, k, ", ")
    end
    print(io, "\n")
end

abstract type AbstractPlanet end
"""
    Planet([derived,] priors, [astrometry,], name=:symbol)

A planet (or substellar companion) part of a model.
Must be constructed with a block of priors, and optionally
additional derived parameters and/or astrometry.
`name` must be a symbol, e.g. `:b`.
"""
struct Planet{TD<:Union{Derived,Nothing},TP<:Union{Priors,Nothing},TA<:Union{Astrometry,Nothing},TPhot<:Union{Photometry,Nothing}} <: AbstractPlanet
    deterministic::TD
    priors::TP
    astrometry::TA
    photometry::TPhot
    name::Symbol
end
export Planet
# There has got to be a better way...
Planet(det::Derived,priors::Priors,astrometry::Union{Astrometry,Nothing}=nothing, photometry::Union{Photometry,Nothing}=nothing; name) = Planet(det,priors, astrometry, photometry, name)
Planet(priors::Priors,astrometry::Union{Astrometry,Nothing}=nothing, photometry::Union{Photometry,Nothing}=nothing; name) = Planet(nothing,priors, astrometry,photometry, name)
Planet(det::Derived,priors::Priors, photometry::Photometry, astrometry::Union{Astrometry,Nothing}=nothing; name) = Planet(det,priors, astrometry, photometry, name)
Planet(priors::Priors, photometry::Photometry, astrometry::Union{Astrometry,Nothing}=nothing; name) = Planet(nothing,priors, astrometry,photometry, name)
Planet(priors::Priors, det::Derived, args...; name) = Planet(det, priors, args...; name)

astrometry(planet::Planet) = planet.astrometry

"""
    System([derived,] priors, [images,] [propermotionanom,] planets..., name=:symbol)

Construct a model of a system.
Must be constructed with a block of priors, and optionally
additional derived parameters.
You may provide `ProperMotionAnom()` and/or `Images()` of the system.
Finally, planet models are listed last.
`name` must be a symbol e.g. `:HD56441`.
"""
struct System{TDet<:Union{Derived,Nothing}, TPriors<:Priors, TPMA<:Union{ProperMotionAnom,Nothing}, TImages<:Union{Nothing,Images},TPlanet}
    deterministic::TDet
    priors::TPriors
    propermotionanom::TPMA
    images::TImages
    planets::TPlanet
    name::Symbol
    function System(
        system_det::Union{Derived,Nothing},
        system_priors::Union{Priors,Nothing},
        propermotionanom::Union{ProperMotionAnom,Nothing},
        images::Union{Images,Nothing},
        planets::AbstractPlanet...;
        name
    )
        if isempty(planets)
            planets_nt = (;)
        else
            planets_nt = namedtuple(
                getproperty.(planets, :name),
                planets
            )
        end
        return new{typeof(system_det), typeof(system_priors), typeof(propermotionanom), typeof(images), typeof(planets_nt)}(
            system_det, system_priors, propermotionanom, images, planets_nt, name
        )
    end
end
export System

# Argument standardization / method cascade.
# Allows users to pass arguments to System in any convenient order.
System(planets::Planet...; kwargs...) = System(Priors(), planets...; kwargs...)
System(priors::Priors, args...; kwargs...) =
    System(nothing, priors, args...,; kwargs...)
System(priors::Priors, planets::Planet...; kwargs...) =
    System(nothing, priors, nothing, nothing, planets...,; kwargs...)
System(priors::Priors, det::Derived, args...; kwargs...) = System(det, priors, args...; kwargs...)
System(det::Derived, priors::Priors, planets::Planet...; kwargs...) =
    System(det, priors, nothing, nothing, planets...; kwargs...)
System(det::Union{Derived,Nothing}, priors::Union{Priors,Nothing}, propermotionanom::ProperMotionAnom, planets::Planet...; kwargs...) =
    System(det, priors, propermotionanom, nothing, planets...; kwargs...)
System(det::Union{Derived,Nothing}, priors::Union{Priors,Nothing}, images::Images, planets::Planet...; kwargs...) =
    System(det, priors, nothing, images, planets...; kwargs...)
System(det::Union{Derived,Nothing}, priors::Union{Priors,Nothing}, images::Images, propermotionanom::ProperMotionAnom, planets::Planet...; kwargs...) =
    System(det, priors, propermotionanom, images, planets...; kwargs...)

#### Show methods
Base.show(io::IO, ::MIME"text/plain", @nospecialize phot::Photometry) = print(
    io, """
        Astrometry[$(length(phot.band))]
        band\tphot\tσ_phot
        ──────────────────────
        $(join(["$(phot.band[i])\t$(phot.phot[i])\t$(phot.σ_phot[i])"
             for i in eachindex(phot.band)],"\n"))
        ──────────────────────
        """)
Base.show(io::IO, ::MIME"text/plain", @nospecialize astrom::Astrometry) = print(
    io, """
        Astrometry[$(length(astrom.epoch))]
        epoch   \tra\tdec\tσ_ra\tσ_dec
        ────────────────────────────────────────────────
        $(join(["$(round(astrom.epoch[i],digits=2))\t$(round(astrom.ra[i],digits=2))\t$(round(astrom.dec[i],digits=2))\t$(round(astrom.σ_ra[i],digits=2))\t$(round(astrom.σ_dec[i],digits=2))\t"
             for i in eachindex(astrom.epoch)],"\n"))
        ────────────────────────────────────────────────
        """)
Base.show(io::IO, ::MIME"text/plain", @nospecialize pma::ProperMotionAnom) = print(
    io, """
        ProperMotionAnom[$(length(pma.ra_epoch))]
        ra_epoch   \tdec_epoch   \tpm_ra\tpm_dec\tσ_pm_ra\tσ_pm_dec
        ──────────────────────────────────────────────────────────────────
        $(join(["$(round(pma.ra_epoch[i],digits=2)) \t$(round(pma.dec_epoch[i],digits=2))   \t$(round(pma.pm_ra[i],digits=2))\t$(round(pma.pm_dec[i],digits=2))\t$(round(pma.σ_pm_ra[i],digits=2))\t$(round(pma.σ_pm_dec[i],digits=2))\t"
             for i in eachindex(pma.ra_epoch)],"\n"))
        ──────────────────────────────────────────────────────────────────
        """)
Base.show(io::IO, ::MIME"text/plain", @nospecialize is::Images) = print(
    io, """
        Images[$(length(is.image))]
        epoch\tband\tplatescale
        ───────────────────────────
        $(join(["$(round(is.epoch[i],digits=2))\t$(is.band[i])\t$(round(is.platescale[i],digits=2))" for i in eachindex(is.epoch)],"\n"))
        ───────────────────────────
        """)
function Base.show(io::IO, mime::MIME"text/plain", @nospecialize p::AbstractPlanet)
    print(io, "Planet $(p.name)\n")
    # if !isnothing(p.deterministic)
    #     show(io, mime, p.deterministic)
    # end
    # show(io, mime, p.priors)
    # if hasproperty(p, :astrometry) && !isnothing(p.astrometry)
    #     show(io, mime, p.astrometry)
    # end
    # print(io, "\n")
    # println(io)
end
function Base.show(io::IO, mime::MIME"text/plain", @nospecialize sys::System)
    print(io, "System model $(sys.name)\n")
    # show(io, mime, sys.priors)
    # for planet in sys.planets
    #     show(io, mime, planet)
    # end
    # if !isnothing(sys.propermotionanom)
    #     show(io, mime, sys.propermotionanom)
    # end
    # if !isnothing(sys.images)
    #     show(io, mime, sys.images)
    # end
    # println(io)
end

"""
    _list_priors(system::System)

Given a System, return a flat vector of prior distributions for all parameters in the System and
any Planets, according to the same order as `make_arr2nt`.
"""
function _list_priors(system::System)

    priors_vec = []
    # System priors
    for prior_distribution in values(system.priors.priors)
        push!(priors_vec,prior_distribution)
    end

    # Planet priors
    for planet in system.planets
        # for prior_distribution in values(planet.priors.priors)
        for (key, prior_distribution) in zip(keys(planet.priors.priors), values(planet.priors.priors))
            push!(priors_vec, prior_distribution)
        end
    end

    # narrow the type
    return map(identity, priors_vec)
end

"""
    make_arr2nt(system::System)

Returns a function that maps an array of parameter values (e.g. sampled from priors,
sampled from posterior, dual numbers, etc.) to a nested named tuple structure suitable
for passing into `ln_like`.
In the process, this evaluates all Deterministic variables defined in the model.

Example:
```julia
julia> arr2nt = DirectDetections.make_arr2nt(system);
julia> params = DirectDetections.sample_priors(system)
9-element Vector{Float64}:
  1.581678216418196
 29.019567489624023
  1.805551918844831
  4.527607604709967
 -2.4432825071288837
 10.165695048003027
 -2.8667829383306063
  0.0006589349005787781
 -1.2600092725864576
julia> arr2nt(params)
(μ = 1.581678216418196, plx = 29.019567489624023, planets = (B = (τ = 1.805551918844831, ω = 4.527607604709967, i = -2.4432825071288837, Ω = 10.165695048003027, loge = -2.8667829383306063, loga = 0.0006589349005787781, logm = -1.2600092725864576, e = 0.001358992505357737, a = 1.0015184052910453, mass = 0.054952914078000334),))
```
"""
function make_arr2nt(system::System)

    # Roll flat vector of variables e.g. drawn from priors or posterior
    # into a nested NamedTuple structure for e.g. ln_like
    i = 0
    body_sys_priors = Expr[]

    # This function uses meta-programming to unroll all the code at compile time.
    # This is a big performance win, since it avoids looping over all the functions
    # etc. In fact, most user defined e.g. Derived callbacks can get inlined right here.

    # Priors
    for key in keys(system.priors.priors)
        i += 1
        ex = :(
            $key = arr[$i]
        )
        push!(body_sys_priors,ex)
    end

    # Deterministic variables for the system
    body_sys_determ = Expr[]
    if !isnothing(system.deterministic)
        for (key,func) in zip(keys(system.deterministic.variables), values(system.deterministic.variables))
            ex = :(
                $key = $func(sys)
            )
            push!(body_sys_determ,ex)
        end
    end

    # Planets: priors & deterministic variables
    body_planets = Expr[]
    for planet in system.planets

        # Priors
        body_planet_priors = Expr[]
        for key in keys(planet.priors.priors)
            i += 1
            ex = :(
                $key = arr[$i]
            )
            push!(body_planet_priors,ex)
        end

        if isnothing(planet.deterministic)
            ex = :(
                $(planet.name) = (;
                    $(body_planet_priors...)
                )
            )
        # Resolve deterministic vars.
        else
            body_planet_determ = Expr[]
            for (key,func) in zip(keys(planet.deterministic.variables), values(planet.deterministic.variables))
                ex = :(
                    $key = $func(sys_res, (;$(body_planet_priors...)))
                )
                push!(body_planet_determ,ex)
            end

            ex = :(
                $(planet.name) = (;
                    $(body_planet_priors...),
                    $(body_planet_determ...)
                )
            )
        end
        push!(body_planets,ex)
    end

    # Here is the function we return.
    # It maps an array of parameters into our nested named tuple structure
    # Note: eval() would normally work fine here, but sometimes we can hit "world age problemms"
    # The RuntimeGeneratedFunctions package avoids these in all cases.
    func = @RuntimeGeneratedFunction(:(function (arr)
        l = $i
        @boundscheck if length(arr) != l
            error("Expected exactly $l elements in array (got $(length(arr)))")
        end
        # Expand system variables from priors
        sys = (;$(body_sys_priors...))
        # Resolve deterministic system variables
        sys_res = (;
            sys...,
            $(body_sys_determ...)
        )
        # Get resolved planets
        pln = (;$(body_planets...))
        # Merge planets into resolved system
        sys_res_pln = (;sys..., planets=pln)
        return sys_res_pln
    end))

    return func
end

ext: jl


text:
push!(LOAD_PATH, "/mnt/g/Python Programme/BALL.jl/src" )
using BALL
using BenchmarkTools

comp = collectSystems(System("../compBench/1htq.pdb"))[2]
dfs = DataFrameSystem(comp)

function dist1()
    firstChain::Union{Chain,Nothing} = comp.first_child_
    resIter::Residue = firstChain.first_child_
    res1 = Residue()
    res2 = Residue()
    for x in 1:60
        if x == 50
            res1 = resIter
        end

        if x ==60
            res2 = resIter
        end
        resIter = resIter.next_
    end
return distance(res1, res2)
end

function dist2()
    return distance(dfs.residues[50], dfs.residues[60])
end

@show @benchmark dist1()
@show @benchmark dist2()


text:
#
# Algebra of queries.
#

import Base:
    convert,
    getindex,
    show,
    >>

#
# Query interface.
#

abstract type AbstractQuery end

"""
    Query(op, args...)

A query is implemented as a pipeline transformation that preserves pipeline
source.  Specifically, a query takes the input pipeline that maps the *source*
to the *input target* and generates a pipeline that maps the *source* to the
*output target*.

Parameter `op` is a function that performs the transformation; `args` are extra
arguments passed to the function.

The query transforms an input pipeline `p` by invoking `op` with the following
arguments:

    op(env::Environment, q::Pipeline, args...)

The result of `op` must be the output pipeline.
"""
struct Query <: AbstractQuery
    op
    args::Vector{Any}

    Query(op; args::Vector{Any}=Any[]) =
        new(op, args)
end

Query(op, args...) =
    Query(op; args=collect(Any, args))

quoteof(F::Query) =
    quoteof(F.op, F.args)

show(io::IO, F::Query) =
    print_expr(io, quoteof(F))

#
# Navigation sugar.
#

struct Navigation <: AbstractQuery
    __path::Tuple{Vararg{Symbol}}
end

Base.getproperty(nav::Navigation, s::Symbol) =
    let path = getfield(nav, :__path)
        Navigation((path..., s))
    end

show(io::IO, nav::Navigation) =
    let path = getfield(nav, :__path)
        print(io, join((:It, path...), "."))
    end

"""
    It :: AbstractQuery

In a query expression, use `It` to refer to the query's input.

```jldoctest
julia> unitknot[Lift(3) >> (It .+ 1)]
┼───┼
│ 4 │
```

`It` is the identity with respect to query composition.

```jldoctest
julia> unitknot[Lift('a':'c') >> It]
──┼───┼
1 │ a │
2 │ b │
3 │ c │
```

`It` provides a shorthand notation for data navigation using
`Get`, so that `It.a.x` is equivalent to `Get(:a) >> Get(:x)`.

```jldoctest
julia> unitknot[Lift((a=(x=1,y=2),)) >> It.a]
│ a    │
│ x  y │
┼──────┼
│ 1  2 │

julia> unitknot[Lift((a=(x=1,y=2),)) >> It.a.x]
│ x │
┼───┼
│ 1 │
```
"""
const It = Navigation(())

translate(mod::Module, ::Val{:it}) = It

#
# Querying a DataKnot.
#

"""
    db::DataKnot[F::Query; params...] :: DataKnot

Queries `db` with `F`.
"""
getindex(db::DataKnot, F; kws...) =
    query(db, Each(F); kws...)

query(db, F; kws...) =
    query(convert(DataKnot, db), Lift(F), sort(collect(Pair{Symbol,DataKnot}, kws), by=first))

function query(db::DataKnot, F::AbstractQuery, params::Vector{Pair{Symbol,DataKnot}}=Pair{Symbol,DataKnot}[])
    db = pack(db, params)
    q = assemble(shape(db), F)
    db′ = q(db)
    return db′
end

function pack(db::DataKnot, params::Vector{Pair{Symbol,DataKnot}})
    !isempty(params) || return db
    ctx_lbls = first.(params)
    ctx_cols = collect(AbstractVector, cell.(last.(params)))
    ctx_shps = collect(AbstractShape, shape.(last.(params)))
    scp_cell = TupleVector(1, AbstractVector[cell(db), TupleVector(ctx_lbls, 1, ctx_cols)])
    scp_shp = TupleOf(shape(db), TupleOf(ctx_lbls, ctx_shps)) |> IsScope
    return DataKnot(scp_shp, scp_cell)
end

assemble(db::DataKnot, F::AbstractQuery, params::Vector{Pair{Symbol,DataKnot}}=Pair{Symbol,DataKnot}[]) =
    assemble(shape(pack(db, params)), F)

assemble(src::AbstractShape, F::AbstractQuery) =
    optimize(uncover(assemble(nothing, cover(src), F)))

assemble(::Nothing, p::Pipeline, F) =
    assemble(Environment(), p, F)

#
# External syntax.
#
"""
     @query expr

Creates a query object from a specialized path-like notation:

* bare identifiers are translated to navigation with `Get`;
* query combinators, such as `Count(X)`, use lower-case names;
* the period (`.`) is used for query composition (`>>`);
* aggregate queries, such as `Count`, require parentheses;
* records can be constructed using curly brackets, `{}`; and
* functions and operators are lifted automatically.

```jldoctest
julia> @query 2x+1
Lift(+, (Lift(*, (Lift(2), Get(:x))), Lift(1)))
```
"""
macro query(ex)
    return :( translate($__module__, $(Expr(:quote, ex)) ) )
end

"""
     @query dataset expr param=...

Applies the query to a dataset with a given set of parameters.

```jldoctest
julia> @query unitknot 2x+1 x=1
┼───┼
│ 3 │
```
"""
macro query(db, exs...)
    exs = map(exs) do ex
        if Meta.isexpr(ex, :(=), 2)
            esc(Expr(:kw, ex.args...))
        else
            :( Each(translate($__module__, $(Expr(:quote, ex)) )) )
        end
    end
    return quote
        query($(esc(db)), $(exs...))
    end
end

function translate(mod::Module, ex::Expr)::AbstractQuery
    head = ex.head
    args = ex.args
    if head === :block
        return Compose(translate.(Ref(mod), filter(arg -> !(arg isa LineNumberNode), args))...)
    elseif head === :. && length(args) == 2 && Meta.isexpr(args[2], :tuple, 1)
        return Compose(translate(mod, args[1]), translate(mod, args[2].args[1]))
    elseif head === :.
        return Compose(translate.(Ref(mod), args)...)
    elseif head === :let && length(args) == 2
        return Given(translate.(Ref(mod), Meta.isexpr(args[1], :block) ? args[1].args : (args[1],))...,
                     translate(mod, args[2]))
    elseif head === :braces
        return Record(translate.(Ref(mod), args)...)
    elseif head === :quote && length(args) == 1 && Meta.isexpr(args[1], :braces)
        return Record(translate.(Ref(mod), args[1].args)...)
    elseif head === :curly && length(args) >= 1
        return Compose(translate(mod, args[1]), Record(translate.(Ref(mod), args[2:end])...))
    elseif head === :call && length(args) >= 1
        call = args[1]
        if call isa QuoteNode
            call = call.value
        end
        if call === :(=>) && length(args) == 3 && args[2] isa Symbol
            return Compose(translate(mod, args[3]), Label(args[2]))
        elseif call isa Symbol
            return translate(mod, Val(call), (args[2:end]...,))
        elseif Meta.isexpr(call, :.) && !isempty(call.args)
            return Compose(translate.(Ref(mod), call.args[1:end-1])...,
                           translate(mod, Expr(:call, call.args[end], args[2:end]...)))
        elseif call isa Base.Callable
            return Lift(call(translate.(Ref(mod), args[2:end])...))
        end
    elseif head == :comparison && length(args) == 3
        return translate(mod, Expr(:call, args[2], args[1], args[3]))
    elseif head == :comparison && length(args) > 3
        return translate(mod, Expr(:&&, Expr(:call, args[2], args[1], args[3]), Expr(head, args[3:end]...)))
    elseif head == :&&
        return Lift(&, (translate.(Ref(mod), args)...,))
    elseif head == :||
        return Lift(|, (translate.(Ref(mod), args)...,))
    end
    error("invalid query expression: $(repr(ex))")
end

translate(mod::Module, sym::Symbol) =
    translate(mod, Val(sym))

translate(::Module, ::Val{:nothing}) =
    Lift(nothing)

translate(::Module, ::Val{:missing}) =
    Lift(missing)

translate(mod::Module, qn::QuoteNode) =
    translate(mod, qn.value)

translate(mod::Module, @nospecialize(v::Val{N})) where {N} =
    Get(N)

function translate(mod::Module, @nospecialize(v::Val{N}), args::Tuple) where {N}
    fn = getfield(mod, N)
    Lift(fn, translate.(Ref(mod), args))
end

translate(mod::Module, val) =
    Lift(val)

#
# Compiling a query.
#

"""
    Environment()

Query compilation state.
"""
mutable struct Environment
end

assemble(env::Environment, p::Pipeline, F)::Pipeline =
    assemble(env, p, Lift(F))

assemble(env::Environment, p::Pipeline, F::Query)::Pipeline =
    F.op(env, p, F.args...)

function assemble(env::Environment, p::Pipeline, nav::Navigation)::Pipeline
    for name in getfield(nav, :__path)
        p = Get(env, p, name)
    end
    p
end

#
# Adapters.
#

# The underlying data shape below flow and scope containers.

domain(shp::AbstractShape) =
    shp

domain(shp::IsFlow) =
    domain(elements(shp))

domain(shp::IsScope) =
    domain(column(shp))

replace_domain(shp::AbstractShape, f) =
    f isa AbstractShape ? f : f(shp)

replace_domain(shp::IsFlow, f) =
    replace_elements(shp, elts -> replace_domain(elts, f))

replace_domain(shp::IsScope, f) =
    replace_column(shp, col -> replace_domain(col, f))

# Finds the output label.

getlabel(p::Pipeline, default) =
    getlabel(target(p), default)

getlabel(shp::AbstractShape, default) =
    default

getlabel(shp::IsLabeled, default) =
    label(shp)

getlabel(shp::IsFlow, default) =
    getlabel(elements(shp), default)

getlabel(shp::IsScope, default) =
    getlabel(column(shp), default)

# Reassigns the output label.

relabel(p::Pipeline, lbl::Union{Symbol,Nothing}) =
    p |> designate(source(p), relabel(target(p), lbl))

relabel(shp::AbstractShape, ::Nothing) =
    shp

relabel(shp::AbstractShape, lbl::Symbol) =
    shp |> IsLabeled(lbl)

relabel(shp::IsLabeled, ::Nothing) =
    subject(shp)

relabel(shp::IsLabeled, lbl::Symbol) =
    subject(shp) |> IsLabeled(lbl)

relabel(shp::IsFlow, lbl::Symbol) =
    replace_elements(shp, elts -> relabel(elts, lbl))

relabel(shp::IsFlow, ::Nothing) =
    replace_elements(shp, elts -> relabel(elts, nothing))

relabel(shp::IsScope, lbl::Symbol) =
    replace_column(shp, col -> relabel(col, lbl))

relabel(shp::IsScope, ::Nothing) =
    replace_column(shp, col -> relabel(col, nothing))

# Removes the flow annotation and strips the scope container from the query output.

function uncover(p::Pipeline)
    q = uncover(target(p))
    chain_of(p, q) |> designate(source(p), target(q))
end

function uncover(src::IsFlow)
    p = uncover(elements(src))
    lbl = getlabel(p, nothing)
    p = relabel(p, nothing)
    tgt = relabel(replace_domain(target(p), dom -> BlockOf(dom, cardinality(src))), lbl)
    with_elements(p) |> designate(src, tgt)
end

uncover(src::IsScope) =
    column(1) |> designate(src, column(src))

uncover(src::AbstractShape) =
    pass() |> designate(src, src)

# Finds or creates a flow container and clones the scope container.

function cover(cell::BlockVector, sig::Signature)
    elts = elements(cell)
    card = cardinality(cell)
    q = elts isa Vector{Union{}} && card == x0to1 ? null_filler() : block_filler(elts, card)
    cover(q |> designate(sig))
end

cover(cell::AbstractVector, sig::Signature) =
    cover(filler(cell[1]) |> designate(sig))

cover(p::Pipeline) =
    cover(source(p), p)

function cover(src::IsScope, p::Pipeline)
    ctx = context(src)
    tgt = TupleOf(target(p), ctx) |> IsScope
    p = tuple_of(p, column(2)) |> designate(src, tgt)
    cover(nothing, p)
end

cover(::AbstractShape, p::Pipeline) =
    cover(nothing, p)

function cover(::Nothing, p::Pipeline)
    q = cover(target(p))
    chain_of(p, q) |> designate(source(p), target(q))
end

cover(src::AbstractShape) =
    wrap() |> designate(src, BlockOf(src, x1to1) |> IsFlow)

cover(src::BlockOf) =
    pass() |> designate(src, src |> IsFlow)

function cover(src::ValueOf)
    ty = eltype(src)
    if ty <: AbstractVector
        ty′ = eltype(ty)
        adapt_vector() |> designate(src, BlockOf(ty′, x0toN) |> IsFlow)
    elseif Missing <: ty
        ty′ = Base.nonmissingtype(ty)
        adapt_missing() |> designate(src, BlockOf(ty′, x0to1) |> IsFlow)
    else
        wrap() |> designate(src, BlockOf(src, x1to1) |> IsFlow)
    end
end

function cover(src::IsLabeled)
    p = cover(subject(src))
    tgt = replace_elements(target(p), IsLabeled(label(src)))
    p |> designate(src, tgt)
end

cover(src::IsFlow) =
    pass() |> designate(src, src)

function cover(src::IsScope)
    p = cover(column(src))
    tgt = target(p)
    tgt = replace_elements(tgt, TupleOf(elements(tgt), context(src)) |> IsScope)
    chain_of(with_column(1, p), distribute(1)) |> designate(src, tgt)
end

#
# Elementwise composition.
#

# Trivial pipes at the source and target endpoints of a pipeline.

source_pipe(p::Pipeline) =
    trivial_pipe(source(p))

target_pipe(p::Pipeline) =
    trivial_pipe(target(p))

trivial_pipe(src::IsFlow) =
    trivial_pipe(elements(src))

trivial_pipe(src::AbstractShape) =
    cover(src)

trivial_pipe(db::DataKnot) =
    cover(shape(db))

# Align pipelines for composition.

realign(p::Pipeline, ::AbstractShape) =
    p

realign(p::Pipeline, ref::IsScope) =
    realign(p, target(p), ref)

realign(p::Pipeline, ::AbstractShape, ::IsScope) =
    p

realign(p::Pipeline, tgt::IsFlow, ref::IsScope) =
    realign(p, elements(tgt), tgt, ref)

realign(p::Pipeline, ::IsScope, ::IsFlow, ::IsScope) =
    p

function realign(p::Pipeline, elts::AbstractShape, tgt::IsFlow, ref::IsScope)
    p′ = chain_of(with_column(1, p), distribute(1))
    ctx = context(ref)
    src′ = TupleOf(source(p), ctx) |> IsScope
    tgt′ = replace_elements(tgt, elts -> TupleOf(elts, ctx) |> IsScope)
    p′ |> designate(src′, tgt′)
end

realign(::AbstractShape, p::Pipeline) =
    p

realign(ref::IsFlow, p::Pipeline) =
    realign(ref, source(p), p)

realign(::IsFlow, ::IsFlow, p::Pipeline) =
    p

realign(ref::IsFlow, ::AbstractShape, p::Pipeline) =
    realign(ref, p, target(p))

function realign(ref::IsFlow, p::Pipeline, ::AbstractShape)
    p′ = with_elements(p)
    src′ = replace_elements(ref, source(p))
    tgt′ = replace_elements(ref, target(p))
    p′ |> designate(src′, tgt′)
end

function realign(ref::IsFlow, p::Pipeline, tgt::IsFlow)
    p′ = chain_of(with_elements(p), flatten())
    src′ = replace_elements(ref, source(p))
    card′ = cardinality(ref)|cardinality(tgt)
    tgt′ = BlockOf(elements(tgt), card′) |> IsFlow
    p′ |> designate(src′, tgt′)
end

# Composition.

compose(p::Pipeline) = p

compose(p1::Pipeline, p2::Pipeline, p3::Pipeline, ps::Pipeline...) =
    foldl(compose, ps, init=compose(compose(p1, p2), p3))

function compose(p1::Pipeline, p2::Pipeline)
    p1 = realign(p1, source(p2))
    p2 = realign(target(p1), p2)
    @assert fits(target(p1), source(p2)) "cannot fit\n$(target(p1))\ninto\n$(source(p2))"
    chain_of(p1, p2) |> designate(source(p1), target(p2))
end

>>(X::Union{DataKnot,AbstractQuery,Pair{Symbol,<:Union{DataKnot,AbstractQuery}}}, Xs...) =
    Compose(X, Xs...)

Compose(Xs...) =
    Query(Compose, Xs...)

quoteof(::typeof(Compose), args::Vector{Any}) =
    quoteof(>>, args)

function Compose(env::Environment, p::Pipeline, Xs...)
    for X in Xs
        p = assemble(env, p, X)
    end
    p
end

#
# Record combinator.
#

function assemble_record(p::Pipeline, xs::Vector{Pipeline})
    lbls = Symbol[]
    cols = Pipeline[]
    seen = Dict{Symbol,Int}()
    for (i, x) in enumerate(xs)
        x = uncover(x)
        lbl = getlabel(x, nothing)
        if lbl !== nothing
            x = relabel(x, nothing)
        else
            lbl = ordinal_label(i)
        end
        if lbl in keys(seen)
            lbls[seen[lbl]] = ordinal_label(seen[lbl])
        end
        seen[lbl] = i
        push!(lbls, lbl)
        push!(cols, x)
    end
    src = elements(target(p))
    tgt = TupleOf(lbls, target.(cols))
    lbl = getlabel(p, nothing)
    if lbl !== nothing
        tgt = relabel(tgt, lbl)
    end
    q = tuple_of(lbls, cols) |> designate(src, tgt)
    q = cover(q)
    compose(p, q)
end

"""
    Record(X₁, X₂ … Xₙ) :: Query

This query emits a record with fields generated by `X₁`, `X₂` … `Xₙ`.

```jldoctest
julia> unitknot[Lift(1:3) >> Record(It, It .* It)]
  │ #A  #B │
──┼────────┼
1 │  1   1 │
2 │  2   4 │
3 │  3   9 │
```

Field labels are inherited from queries.

```jldoctest
julia> unitknot[Lift(1:3) >> Record(:x => It,
                                    :x² => It .* It)]
  │ x  x² │
──┼───────┼
1 │ 1   1 │
2 │ 2   4 │
3 │ 3   9 │
```
"""
Record(Xs...) =
    Query(Record, Xs...)

function Record(env::Environment, p::Pipeline, Xs...)
    xs = assemble.(Ref(env), Ref(target_pipe(p)), collect(AbstractQuery, Xs))
    assemble_record(p, xs)
end

translate(mod::Module, ::Val{:record}, args::Tuple) =
    Record(translate.(Ref(mod), args)...)

#
# Mix combinator.
#

function assemble_mix(p::Pipeline, xs::Vector{Pipeline})
    lbls = Symbol[]
    cols = Pipeline[]
    card = x1to1
    seen = Dict{Symbol,Int}()
    for (i, x) in enumerate(xs)
        card |= cardinality(target(x))
        x = uncover(x)
        lbl = getlabel(x, nothing)
        if lbl !== nothing
            x = relabel(x, nothing)
        else
            lbl = ordinal_label(i)
        end
        if lbl in keys(seen)
            lbls[seen[lbl]] = ordinal_label(seen[lbl])
        end
        seen[lbl] = i
        push!(lbls, lbl)
        x = chain_of(
                x,
                with_elements(wrap()),
        ) |> designate(source(x), replace_elements(target(x), elts -> BlockOf(elts, x1to1)))
        push!(cols, x)
    end
    src = elements(target(p))
    tgt = BlockOf(TupleOf(lbls, elements.(target.(cols))), card)
    lbl = getlabel(p, nothing)
    if lbl !== nothing
        tgt = relabel(tgt, lbl)
    end
    q = chain_of(
            tuple_of(lbls, cols),
            distribute_all(),
    ) |> designate(src, tgt)
    q = cover(q)
    compose(p, q)
end

"""
    Mix(X₁, X₂ … Xₙ) :: Query

This query emits records containing every combination of elements
generated by `X₁`, `X₂` … `Xₙ`.

```jldoctest
julia> unitknot[Mix(Lift(1:2), Lift('a':'c'))]
  │ #A  #B │
──┼────────┼
1 │  1  a  │
2 │  1  b  │
3 │  1  c  │
4 │  2  a  │
5 │  2  b  │
6 │  2  c  │
```
"""
Mix(Xs...) =
    Query(Mix, Xs...)

function Mix(env::Environment, p::Pipeline, Xs...)
    xs = assemble.(Ref(env), Ref(target_pipe(p)), collect(AbstractQuery, Xs))
    assemble_mix(p, xs)
end

translate(mod::Module, ::Val{:mix}, args::Tuple) =
    Mix(translate.(Ref(mod), args)...)

#
# Collect combinator.
#

function as_record(p::Pipeline)
    q = as_record(elements(target(p)))
    q !== nothing ?
        compose(p, cover(q)) :
        p
end

as_record(src::AbstractShape) = nothing

as_record(::TupleOf) = nothing

function as_record(src::IsLabeled)
    p = as_record(subject(src))
    p !== nothing ?
        p |> designate(src, target(p) |> IsLabeled(label(src))) :
        nothing
end

as_record(src::IsFlow) =
    as_record(elements(src))

function as_record(src::IsScope)
    p = as_record(column(src))
    p !== nothing ?
        chain_of(column(1), p) |> designate(src, target(p)) :
        nothing
end

as_record(src::ValueOf) =
    as_record(src.ty)

as_record(::Type) =
    nothing

function as_record(ity::Type{<:NamedTuple})
    lbls = collect(Symbol, ity.parameters[1])
    cols = collect(AbstractShape, ity.parameters[2].parameters)
    adapt_tuple() |> designate(ity, TupleOf(lbls, cols))
end

function as_record(ity::Type{<:Tuple})
    cols = collect(AbstractShape, ity.parameters)
    adapt_tuple() |> designate(ity, TupleOf(cols))
end

function assemble_collect(p, x)
    p = as_record(p)
    src = elements(target(p))
    dom = deannotate(domain(src))
    dom isa TupleOf || error("expected a record; got\n$(syntaxof(dom))")
    x = uncover(x)
    x_lbl = getlabel(x, nothing)
    x = relabel(x, nothing)
    cols = Pipeline[]
    lbls = Symbol[]
    x_pos = width(dom)+1
    for i in 1:width(dom)
        lbl = label(dom, i)
        if lbl == x_lbl
            x_pos = i
            continue
        end
        col = lookup(src, i)
        push!(cols, col)
        if lbl == ordinal_label(i)
            lbl = ordinal_label(length(cols))
        end
        push!(lbls, lbl)
    end
    if !fits(target(x), BlockOf(Nothing, x1to1))
        splice!(cols, x_pos:x_pos-1, Ref(x))
        splice!(lbls, x_pos:x_pos-1, Ref(x_lbl !== nothing ? x_lbl : ordinal_label(length(cols))))
    end
    tgt = TupleOf(lbls, target.(cols))
    lbl = getlabel(p, nothing)
    if lbl !== nothing
        tgt = relabel(tgt, lbl)
    end
    q = tuple_of(lbls, cols) |> designate(src, tgt)
    q = cover(q)
    compose(p, q)
end

"""
    Collect(X₁, X₂ … Xₙ) :: Query

In the combinator form, `Collect(X₁, X₂ … Xₙ)` adds fields `X₁`,
`X₂` … `Xₙ` to the input record.

```jldoctest
julia> unitknot[Record(:x => 1) >> Collect(:y => 2 .* It.x)]
│ x  y │
┼──────┼
│ 1  2 │
```

If a field already exists, it is replaced.

```jldoctest
julia> unitknot[Record(:x => 1) >> Collect(:x => 2 .* It.x)]
│ x │
┼───┼
│ 2 │
```

To remove a field, assign it the value `nothing`.

```jldoctest
julia> unitknot[Record(:x => 1) >> Collect(:y => 2 .* It.x, :x => nothing)]
│ y │
┼───┼
│ 2 │
```

---

    Each(X >> Record) :: Query

In the query form, `Collect` appends a field to the source record.

```jldoctest
julia> unitknot[Lift(1) >> Label(:x) >> Collect]
│ x │
┼───┼
│ 1 │
```
"""
Collect(X, Ys...) =
    Query(Collect, X, Ys...)

Lift(::typeof(Collect)) =
    Then(Collect)

Collect(env::Environment, p::Pipeline, X, Ys...) =
    Collect(env, Collect(env, p, X), Ys...)

function Collect(env::Environment, p::Pipeline, X)
    x = assemble(env, target_pipe(p), X)
    assemble_collect(p, x)
end

translate(mod::Module, ::Val{:collect}, args::Tuple) =
    Collect(translate.(Ref(mod), args)...)

translate(mod::Module, ::Val{:collect}, ::Tuple{}) =
    Then(Collect)

#
# Join combinator.
#

join_pipe(p::Pipeline) =
    join_pipe(source(p), target(p))

join_pipe(src::AbstractShape, dst::AbstractShape) =
    trivial_pipe(replace_domain(dst, domain(src)))

function assemble_join(p::Pipeline, x::Pipeline)
    p = as_record(p)
    dom = deannotate(domain(target(p)))
    dom isa TupleOf || error("expected a record; got\n$(syntaxof(dom))")
    p0 = uncover(source(p))
    q = assemble_join(target(p), target(p0), x)
    chain_of(tuple_of(p, p0), q) |> designate(source(p), target(q))
end

function assemble_join(src::IsFlow, src0::AbstractShape, x::Pipeline)
    q = assemble_join(elements(src), src0, x)
    q′ = chain_of(distribute(1), with_elements(q))
    q′ |> designate(TupleOf(src, src0), replace_elements(src, target(q)))
end

function assemble_join(src::IsScope, src0::AbstractShape, x::Pipeline)
    q = assemble_join(column(src), replace_column(src, src0), x)
    q′ = tuple_of(chain_of(tuple_of(chain_of(column(1), column(1)),
                                    tuple_of(column(2),
                                             chain_of(column(1), column(2)))),
                           q),
                  chain_of(chain_of(column(1), column(2))))
    q′ |> designate(TupleOf(src, src0), replace_column(src, target(q)))
end

function assemble_join(src::AbstractShape, src0::AbstractShape, x::Pipeline)
    dom = deannotate(domain(src))::TupleOf
    x = uncover(x)
    x_lbl = getlabel(x, nothing)
    x = relabel(x, nothing)
    cols = Pipeline[]
    col_shps = AbstractShape[]
    lbls = Symbol[]
    x_pos = width(dom)+1
    for i in 1:width(dom)
        lbl = label(dom, i)
        if lbl == x_lbl
            x_pos = i
            continue
        end
        col = lookup(src, i)
        push!(cols, chain_of(column(1), col))
        push!(col_shps, target(col))
        if lbl == ordinal_label(i)
            lbl = ordinal_label(length(cols))
        end
        push!(lbls, lbl)
    end
    splice!(cols, x_pos:x_pos-1, Ref(chain_of(column(2), x)))
    splice!(col_shps, x_pos:x_pos-1, Ref(target(x)))
    splice!(lbls, x_pos:x_pos-1, Ref(x_lbl !== nothing ? x_lbl : ordinal_label(length(cols))))
    tgt = TupleOf(lbls, col_shps)
    lbl = getlabel(src, nothing)
    if lbl !== nothing
        tgt = relabel(tgt, lbl)
    end
    tuple_of(lbls, cols) |> designate(TupleOf(src, src0), tgt)
end

"""
    Join(X) :: Query

`Join(X)` evaluates `X` in the source context and adds it
as a field to the input record.

```jldoctest
julia> unitknot[Record(:x => 1) >> Each(Record(:y => 2 .* It.x) >> Join(It.x))]
│ y  x │
┼──────┼
│ 2  1 │
```
"""
Join(X) =
    Query(Join, X)

function Join(env::Environment, p::Pipeline, X)
    x = assemble(env, join_pipe(p), X)
    assemble_join(p, x)
end

translate(mod::Module, ::Val{:join}, (arg,)::Tuple{Any}) =
    Join(translate(mod, arg))

#
# Lifting Julia values and functions.
#

function assemble_lift(p::Pipeline, f, xs::Vector{Pipeline})
    cols = uncover.(xs)
    ity = Tuple{eltype.(target.(cols))...}
    oty = Core.Compiler.return_type(f, ity)
    oty != Union{} || error("cannot apply $f to $ity")
    src = elements(target(p))
    tgt = ValueOf(oty)
    q = if length(cols) == 1
            card = cardinality(target(xs[1]))
            if fits(x1toN, card) && !(oty <: AbstractVector)
                chain_of(cols[1], block_lift(f))
            else
                chain_of(cols[1], lift(f))
            end
        else
            chain_of(tuple_of(Symbol[], cols), tuple_lift(f))
        end |> designate(src, tgt)
    q = cover(q)
    compose(p, q)
end

Lift(X::AbstractQuery) = X

"""
    Lift(val) :: Query

This converts any value to a constant query.

```jldoctest
julia> unitknot[Lift("Hello")]
┼───────┼
│ Hello │
```

`AbstractVector` objects become plural queries.

```jldoctest
julia> unitknot[Lift('a':'c')]
──┼───┼
1 │ a │
2 │ b │
3 │ c │
```

To specify the vector cardinality, add `:x0to1`, `:x0toN`,
`:x1to1`, or `:x1toN`.

```jldoctest
julia> unitknot[Lift('a':'c', :x1toN)]
──┼───┼
1 │ a │
2 │ b │
3 │ c │
```

The `missing` value makes an query with no output.

```jldoctest
julia> unitknot[Lift(missing)]
(empty)
```
"""
Lift(val) =
    Query(Lift, val)

Lift(elts::AbstractVector, card::Union{Cardinality,Symbol}) =
    Query(Lift, elts, card)

"""
    Lift(f, (X₁, X₂ … Xₙ)) :: Query

`Lift` lets you use a function as a query combinator.

```jldoctest
julia> unitknot[Lift((x=1, y=2)) >> Lift(+, (It.x, It.y))]
┼───┼
│ 3 │
```

`Lift` is implicitly used when a function is broadcast over
queries.

```jldoctest
julia> unitknot[Lift((x=1, y=2)) >> (It.x .+ It.y)]
┼───┼
│ 3 │
```

Functions accepting a `AbstractVector` can be used with plural
queries.

```jldoctest
julia> unitknot[sum.(Lift(1:3))]
┼───┼
│ 6 │
```

Functions returning `AbstractVector` become plural queries.

```jldoctest
julia> unitknot[Lift((x='a', y='c')) >> Lift(:, (It.x, It.y))]
──┼───┼
1 │ a │
2 │ b │
3 │ c │
```
"""
Lift(f, Xs::Tuple) =
    Query(Lift, f, Xs)

convert(::Type{AbstractQuery}, val) =
    Lift(val)

convert(::Type{AbstractQuery}, F::AbstractQuery) =
    F

Lift(env::Environment, p::Pipeline, val) =
    Lift(env, p, convert(DataKnot, val))

Lift(env::Environment, p::Pipeline, elts::AbstractVector, card::Union{Cardinality,Symbol}) =
    Lift(env, p, DataKnot(Any, elts, card))

function Lift(env::Environment, p::Pipeline, f, Xs::Tuple)
    xs = assemble.(Ref(env), Ref(target_pipe(p)), collect(AbstractQuery, Xs))
    assemble_lift(p, f, xs)
end

function Lift(env::Environment, p::Pipeline, db::DataKnot)
    q = cover(cell(db), Signature(elements(target(p)), shape(db)))
    compose(p, q)
end

# Broadcasting.

struct QueryStyle <: Base.BroadcastStyle
end

Base.BroadcastStyle(::Type{<:Union{AbstractQuery,DataKnot,Pair{Symbol,<:Union{AbstractQuery,DataKnot}}}}) =
    QueryStyle()

Base.BroadcastStyle(s::QueryStyle, ::Broadcast.DefaultArrayStyle) =
    s

Base.broadcastable(X::Union{AbstractQuery,DataKnot,Pair{Symbol,<:Union{AbstractQuery,DataKnot}}}) =
    X

Base.Broadcast.instantiate(bc::Broadcast.Broadcasted{QueryStyle}) =
    bc

Base.copy(bc::Broadcast.Broadcasted{QueryStyle}) =
    BroadcastLift(bc)

BroadcastLift(bc::Broadcast.Broadcasted{QueryStyle}) =
    BroadcastLift(bc.f, (BroadcastLift.(bc.args)...,))

BroadcastLift(val) = val

BroadcastLift(f, Xs) = Query(BroadcastLift, f, Xs)

BroadcastLift(env::Environment, p::Pipeline, args...) =
    Lift(env, p, args...)

quoteof(::typeof(BroadcastLift), args::Vector{Any}) =
    quoteof(broadcast, Any[args[1], quoteof.(args[2])...])

Lift(bc::Broadcast.Broadcasted{QueryStyle}) =
    BroadcastLift(bc)

#
# Each combinator.
#

"""
    Each(X) :: Query

This evaluates `X` elementwise.

```jldoctest
julia> X = Lift('a':'c') >> Count;

julia> unitknot[Lift(1:3) >> Each(X)]
──┼───┼
1 │ 3 │
2 │ 3 │
3 │ 3 │
```

Compare this with the query without `Each`.

```jldoctest
julia> X = Lift('a':'c') >> Count;

julia> unitknot[Lift(1:3) >> X]
┼───┼
│ 9 │
```
"""
Each(X) = Query(Each, X)

Each(env::Environment, p::Pipeline, X) =
    compose(p, assemble(env, target_pipe(p), X))

translate(mod::Module, ::Val{:each}, (arg,)::Tuple{Any}) =
    Each(translate(mod, arg))

#
# Assigning labels.
#

"""
    Label(lbl::Symbol) :: Query

This assigns a label to the output.

```jldoctest
julia> unitknot[Lift("Hello World") >> Label(:greeting)]
│ greeting    │
┼─────────────┼
│ Hello World │
```

A label could also be assigned using the `=>` operator.

```jldoctest
julia> unitknot[:greeting => Lift("Hello World")]
│ greeting    │
┼─────────────┼
│ Hello World │
```
"""
Label(lbl::Symbol) =
    Query(Label, lbl)

Label(env::Environment, p::Pipeline, lbl::Symbol) =
    relabel(p, lbl)

Lift(p::Pair{Symbol}) =
    Compose(p.second, Label(p.first))

translate(mod::Module, ::Val{:label}, (arg,)::Tuple{Symbol}) =
    Label(arg)

#
# Assigning a name to a query.
#

"""
    Tag(name::Symbol, F) :: Query

This provides a substitute name for a query.

```jldoctest
julia> IncIt = It .+ 1
It .+ 1

julia> IncIt = Tag(:IncIt, It .+ 1)
IncIt
```

---

    Tag(name::Symbol, (X₁, X₂ … Xₙ), F) :: Query

This provides a substitute name for a query combinator.

```jldoctest
julia> Inc(X) = Lift(+, (X, 1));

julia> Inc(It)
Lift(+, (It, 1))

julia> Inc(X) = Tag(:Inc, (X,), Lift(+, (X, 1)));

julia> Inc(It)
Inc(It)
```
"""
Tag(name::Symbol, X) =
    Query(Tag, name, X)

Tag(name::Symbol, args::Tuple, X) =
    Query(Tag, name, args, X)

Tag(F::Union{Function,DataType}, args::Tuple, X) =
    Tag(nameof(F), args, X)

Tag(env::Environment, p::Pipeline, name, X) =
    assemble(env, p, X)

Tag(env::Environment, p::Pipeline, name, args, X) =
    assemble(env, p, X)

quoteof(::typeof(Tag), args::Vector{Any}) =
    quoteof(Tag, args...)

quoteof(::typeof(Tag), name::Symbol, X) =
    name

quoteof(::typeof(Tag), name::Symbol, args::Tuple, X) =
    Expr(:call, name, quoteof.(args)...)

#
# Attributes and parameters.
#

"""
    Get(lbl::Symbol) :: Query

This query extracts a field value by its label.

```jldoctest
julia> unitknot[Lift((x=1, y=2)) >> Get(:x)]
│ x │
┼───┼
│ 1 │
```

This has a shorthand form using `It`.

```jldoctest
julia> unitknot[Lift((x=1, y=2)) >> It.x]
│ x │
┼───┼
│ 1 │
```

With unlabeled fields, ordinal labels (A, B, ...) can be used.

```jldoctest
julia> unitknot[Lift((1,2)) >> It.B]
┼───┼
│ 2 │
```
"""
Get(name) =
    Query(Get, name)

function Get(env::Environment, p::Pipeline, name)
    tgt = target(p)
    q = lookup(tgt, name)
    q !== nothing || error("cannot find \"$name\" at\n$(syntaxof(tgt))")
    q = cover(q)
    compose(p, q)
end

lookup(::AbstractShape, ::Any) = nothing

lookup(src::IsLabeled, name::Any) =
    lookup(subject(src), name)

lookup(src::IsFlow, name::Any) =
    lookup(elements(src), name)

function lookup(src::IsScope, name::Any)
    q = lookup(column(src), name)
    q !== nothing ?
        chain_of(column(1), q) |> designate(src, target(q)) :
        nothing
end

function lookup(src::IsScope, name::Symbol)
    q = lookup(context(src), name)
    q === nothing || return chain_of(column(2), q) |> designate(src, target(q))
    q = lookup(column(src), name)
    q === nothing || return chain_of(column(1), q) |> designate(src, target(q))
    nothing
end

function lookup(lbls::Vector{Symbol}, name::Symbol)
    j = findlast(isequal(name), lbls)
    if j === nothing
        j = findlast(isequal(Symbol("#$name")), lbls)
    end
    j
end

lookup(src::TupleOf, j::Int) =
    column(j) |> designate(src, column(src, j))

function lookup(src::TupleOf, name::Symbol)
    lbls = labels(src)
    if isempty(lbls)
        lbls = Symbol[ordinal_label(i) for i = 1:width(src)]
    end
    j = lookup(lbls, name)
    j !== nothing || return nothing

    tgt = relabel(column(src, j), name == lbls[j] ? name : nothing)
    column(lbls[j]) |> designate(src, tgt)
end

lookup(src::ValueOf, name) =
    lookup(src.ty, name)

lookup(::Type, ::Any) =
    nothing

function lookup(ity::Type{<:NamedTuple}, name::Symbol)
    j = lookup(collect(Symbol, ity.parameters[1]), name)
    j !== nothing || return nothing
    oty = ity.parameters[2].parameters[j]
    lift(getindex, j) |> designate(ity, oty |> IsLabeled(name))
end

function lookup(ity::Type{<:Tuple}, name::Symbol)
    lbls = Symbol[ordinal_label(i) for i = 1:length(ity.parameters)]
    j = lookup(lbls, name)
    j !== nothing || return nothing
    oty = ity.parameters[j]
    lift(getindex, j) |> designate(ity, oty)
end

#
# Specifying context parameters.
#

function assemble_keep(p::Pipeline, q::Pipeline)
    q = uncover(q)
    tgt = target(q)
    name = getlabel(tgt, nothing)
    name !== nothing || error("parameter name is not specified")
    tgt = relabel(tgt, nothing)
    lbls′ = Symbol[]
    cols′ = AbstractShape[]
    perm = Int[]
    src = source(q)
    if src isa IsScope
        ctx = context(src)
        for j = 1:width(ctx)
            lbl = label(ctx, j)
            if lbl != name
                push!(lbls′, lbl)
                push!(cols′, column(ctx, j))
                push!(perm, j)
            end
        end
    end
    push!(lbls′, name)
    push!(cols′, tgt)
    ctx′ = TupleOf(lbls′, cols′)
    qs = Pipeline[chain_of(column(2), column(j)) for j in perm]
    push!(qs, q)
    tgt = BlockOf(TupleOf(src isa IsScope ? column(src) : src, ctx′) |> IsScope,
                  x1to1) |> IsFlow
    q = chain_of(tuple_of(src isa IsScope ? column(1) : pass(),
                          tuple_of(lbls′, qs)),
                 wrap(),
    ) |> designate(src, tgt)
    compose(p, q)
end

"""
    Keep(X₁, X₂ … Xₙ) :: Query

`Keep` evaluates named queries, making their results available for
subsequent computation.

```jldoctest
julia> unitknot[Keep(:x => 2) >> It.x]
│ x │
┼───┼
│ 2 │
```

`Keep` does not otherwise change its input.

```jldoctest
julia> unitknot[Lift(1) >> Keep(:x => 2) >> (It .+ It.x)]
┼───┼
│ 3 │
```
"""
Keep(P, Qs...) =
    Query(Keep, P, Qs...)

Keep(env::Environment, p::Pipeline, P, Qs...) =
    Keep(env, Keep(env, p, P), Qs...)

function Keep(env::Environment, p::Pipeline, P)
    q = assemble(env, target_pipe(p), P)
    assemble_keep(p, q)
end

translate(mod::Module, ::Val{:keep}, args::Tuple{Any,Vararg{Any}}) =
    Keep(translate.(Ref(mod), args)...)

#
# Setting the scope for context parameters.
#

function assemble_given(p::Pipeline, q::Pipeline)
    q = cover(uncover(q))
    compose(p, q)
end

"""
    Given(X₁, X₂ … Xₙ, Q) :: Query

This evaluates `Q` in a context augmented with named parameters
added by a set of queries.

```jldoctest
julia> unitknot[Given(:x => 2, It.x .+ 1)]
┼───┼
│ 3 │
```
"""
Given(P, Xs...) =
    Query(Given, P, Xs...)

Given(env::Environment, p::Pipeline, Xs...) =
    Given(env, p, Keep(Xs[1:end-1]...) >> Each(Xs[end]))

function Given(env::Environment, p::Pipeline, X)
    q = assemble(env, target_pipe(p), X)
    assemble_given(p, q)
end

const Let = Given

translate(mod::Module, ::Val{:given}, args::Tuple{Any,Vararg{Any}}) =
    Given(translate.(Ref(mod), args)...)

#
# Then assembly.
#

Then(q::Pipeline) =
    Query(Then, q)

Then(env::Environment, p::Pipeline, q::Pipeline) =
    compose(p, q)

Then(ctor) =
    Query(Then, ctor)

Then(ctor, args::Tuple) =
    Query(Then, ctor, args)

Then(env::Environment, p::Pipeline, ctor, args::Tuple=()) =
    assemble(env, source_pipe(p), ctor(Then(p), args...))

#
# Count and other aggregate combinators.
#

function assemble_count(p::Pipeline)
    p = uncover(p)
    q = chain_of(p,
                 block_length(),
    ) |> designate(source(p), Int)
    cover(q)
end

function assemble_exists(p::Pipeline)
    p = uncover(p)
    q = chain_of(p,
                 block_not_empty(),
    ) |> designate(source(p), Bool)
    cover(q)
end

"""
    Count(X) :: Query

In the combinator form, `Count(X)` emits the number of elements
produced by `X`.

```jldoctest
julia> X = Lift('a':'c');

julia> unitknot[Count(X)]
┼───┼
│ 3 │
```

---

    Each(X >> Count) :: Query

In the query form, `Count` emits the number of elements in its input.

```jldoctest
julia> X = Lift('a':'c');

julia> unitknot[X >> Count]
┼───┼
│ 3 │
```

To limit the scope of aggregation, use `Each`.

```jldoctest
julia> X = Lift('a':'c');

julia> unitknot[Lift(1:3) >> Each(X >> Count)]
──┼───┼
1 │ 3 │
2 │ 3 │
3 │ 3 │
```
"""
Count(X) =
    Query(Count, X)

Lift(::typeof(Count)) =
    Then(Count)

function Count(env::Environment, p::Pipeline, X)
    x = assemble(env, target_pipe(p), X)
    compose(p, assemble_count(x))
end

"""
    Exists(X) :: Query

In the combinator form, `Exists(X)` emits a boolean testing if `X` produces any elements.

```jldoctest
julia> X = Lift('a':'c');

julia> unitknot[Exists(X)]
┼──────┼
│ true │
```

When the query argument `X` is empty, `Exists(X)` produces `false`.

```jldoctest
julia> X = Lift([]);

julia> unitknot[Exists(X)]
┼───────┼
│ false │
```

---

    Each(X >> Exists) :: Query

In the query form, `Exists` emits a boolean testing if its input has any elements.

```jldoctest
julia> X = Lift('a':'c');

julia> unitknot[X >> Exists]
┼──────┼
│ true │
```

When the query input is empty, `Exists` produces `false`.

```jldoctest
julia> X = Lift([]);

julia> unitknot[X >> Exists]
┼───────┼
│ false │
```
"""
Exists(X) =
    Query(Exists, X)

Lift(::typeof(Exists)) =
    Then(Exists)

function Exists(env::Environment, p::Pipeline, X)
    x = assemble(env, target_pipe(p), X)
    compose(p, assemble_exists(x))
end

"""
    Sum(X) :: Query

In the combinator form, `Sum(X)` emits the sum of elements
produced by `X`.

```jldoctest
julia> X = Lift(1:3);

julia> unitknot[Sum(X)]
┼───┼
│ 6 │
```

The `Sum` of an empty input is `0`.

```jldoctest
julia> unitknot[Sum(Int[])]
┼───┼
│ 0 │
```

---

    Each(X >> Sum) :: Query

In the query form, `Sum` emits the sum of input elements.

```jldoctest
julia> X = Lift(1:3);

julia> unitknot[X >> Sum]
┼───┼
│ 6 │
```
"""
Sum(X) =
    Query(Sum, X)

Lift(::typeof(Sum)) =
    Then(Sum)

"""
     Max(X) :: Query

In the combinator form, `Max(X)` finds the maximum among the
elements produced by `X`.

```jldoctest
julia> X = Lift(1:3);

julia> unitknot[Max(X)]
┼───┼
│ 3 │
```

The `Max` of an empty input is empty.

```jldoctest
julia> unitknot[Max(Int[])]
(empty)
```

---

    Each(X >> Max) :: Query

In the query form, `Max` finds the maximum of its input elements.

```jldoctest
julia> X = Lift(1:3);

julia> unitknot[X >> Max]
┼───┼
│ 3 │
```
"""
Max(X) =
    Query(Max, X)

Lift(::typeof(Max)) =
    Then(Max)

"""
     Min(X) :: Query

In the combinator form, `Min(X)` finds the minimum among the
elements produced by `X`.

```jldoctest
julia> X = Lift(1:3);

julia> unitknot[Min(X)]
┼───┼
│ 1 │
```

The `Min` of an empty input is empty.

```jldoctest
julia> unitknot[Min(Int[])]
(empty)
```

---

    Each(X >> Min) :: Query

In the query form, `Min` finds the minimum of its input elements.

```jldoctest
julia> X = Lift(1:3);

julia> unitknot[X >> Min]
┼───┼
│ 1 │
```
"""
Min(X) =
    Query(Min, X)

Lift(::typeof(Min)) =
    Then(Min)

function Sum(env::Environment, p::Pipeline, X)
    x = assemble(env, target_pipe(p), X)
    assemble_lift(p, sum, Pipeline[x])
end

maximum_missing(v) =
    !isempty(v) ? maximum(v) : missing

function Max(env::Environment, p::Pipeline, X)
    x = assemble(env, target_pipe(p), X)
    card = cardinality(target(x))
    optional = fits(x0to1, card)
    assemble_lift(p, optional ? maximum_missing : maximum, Pipeline[x])
end

minimum_missing(v) =
    !isempty(v) ? minimum(v) : missing

function Min(env::Environment, p::Pipeline, X)
    x = assemble(env, target_pipe(p), X)
    card = cardinality(target(x))
    optional = fits(x0to1, card)
    assemble_lift(p, optional ? minimum_missing : minimum, Pipeline[x])
end

translate(mod::Module, ::Val{:count}, (arg,)::Tuple{Any}) =
    Count(translate(mod, arg))

translate(mod::Module, ::Val{:count}, ::Tuple{}) =
    Then(Count)

translate(mod::Module, ::Val{:exists}, (arg,)::Tuple{Any}) =
    Exists(translate(mod, arg))

translate(mod::Module, ::Val{:exists}, ::Tuple{}) =
    Then(Exists)

translate(mod::Module, ::Val{:sum}, (arg,)::Tuple{Any}) =
    Sum(translate(mod, arg))

translate(mod::Module, ::Val{:sum}, ::Tuple{}) =
    Then(Sum)

translate(mod::Module, ::Val{:max}, (arg,)::Tuple{Any}) =
    Max(translate(mod, arg))

translate(mod::Module, ::Val{:max}, ::Tuple{}) =
    Then(Max)

translate(mod::Module, ::Val{:min}, (arg,)::Tuple{Any}) =
    Min(translate(mod, arg))

translate(mod::Module, ::Val{:min}, ::Tuple{}) =
    Then(Min)

#
# Filter combinator.
#

function assemble_filter(p::Pipeline, x::Pipeline)
    x = uncover(x)
    fits(target(x), BlockOf(ValueOf(Bool))) || error("expected a predicate")
    q = chain_of(tuple_of(pass(),
                          chain_of(x, block_any())),
                 sieve_by(),
    ) |> designate(source(x), BlockOf(source(x), x0to1) |> IsFlow)
    compose(p, q)
end

"""
    Filter(X) :: Query

This query emits the elements from its input that satisfy a given
condition.

```jldoctest
julia> unitknot[Lift(1:5) >> Filter(isodd.(It))]
──┼───┼
1 │ 1 │
2 │ 3 │
3 │ 5 │
```

When the predicate query produces an empty output, the condition
is presumed to have failed.

```jldoctest
julia> unitknot[Lift('a':'c') >> Filter(missing)]
(empty)
```

When the predicate produces plural output, the condition succeeds
if at least one output value is `true`.

```jldoctest
julia> unitknot[Lift('a':'c') >> Filter([true,false])]
──┼───┼
1 │ a │
2 │ b │
3 │ c │
```
"""
Filter(X) =
    Query(Filter, X)

function Filter(env::Environment, p::Pipeline, X)
    x = assemble(env, target_pipe(p), X)
    assemble_filter(p, x)
end

translate(mod::Module, ::Val{:filter}, (arg,)::Tuple{Any}) =
    Filter(translate(mod, arg))

#
# First, Last, Nth.
#

assemble_first(p::Pipeline, inv::Bool=false) =
    assemble_nth(p, !inv ? 1 : -1, cardinality(target(p))&x0to1)

function assemble_nth(p::Pipeline, n::Int, card::Cardinality=x0to1)
    elts = elements(target(p))
    chain_of(
        p,
        get_by(n, card),
    ) |> designate(source(p), BlockOf(elts, card) |> IsFlow)
end

function assemble_nth(p::Pipeline, n::Pipeline)
    n = uncover(n)
    fits(target(n), BlockOf(ValueOf(Int), x1to1)) || error("expected a singular mandatory integer")
    src = source(p)
    tgt = BlockOf(elements(target(p)), x0to1) |> IsFlow
    chain_of(
        tuple_of(p, n),
        get_by(),
    ) |> designate(src, tgt)
end

"""
    First(X) :: Query

In the combinator form, `First(X)` emits the first element produced by its argument `X`.

```jldoctest
julia> X = Lift('a':'c');

julia> unitknot[First(X)]
┼───┼
│ a │
```
---

    Each(X >> First) :: Query

In the query form, `First` emits the first element of its input.

```jldoctest
julia> X = Lift('a':'c');

julia> unitknot[X >> First]
┼───┼
│ a │
```
"""
First(X) = Query(First, X)

"""
    Last(X) :: Query

In the combinator form, `Last(X)` emits the last element produced by its argument `X`.

```jldoctest
julia> X = Lift('a':'c');

julia> unitknot[Last(X)]
┼───┼
│ c │
```
---

    Each(X >> Last) :: Query

In the query form, `Last` emits the last element of its input.

```jldoctest
julia> X = Lift('a':'c');

julia> unitknot[X >> Last]
┼───┼
│ c │
```
"""
Last(X) = Query(Last, X)

"""
    Nth(X, N) :: Query

In the combinator form, `Nth(X, N)` emits the `N`th element produced by its argument `X`.

```jldoctest
julia> X = Lift('a':'d');

julia> N = Count(X) .÷ 2;

julia> unitknot[Nth(X, N)]
┼───┼
│ b │
```
---

    Each(X >> Nth(N)) :: Query

In the query form, `Nth(N)` emits the `N`th element produced by its input.

```jldoctest
julia> X = Lift('a':'d');

julia> N = Count(X) .÷ 2;

julia> unitknot[X >> Nth(N)]
┼───┼
│ b │
```
"""
Nth(X, N) = Query(Nth, X, N)

Lift(::typeof(First)) =
    Then(First)

Lift(::typeof(Last)) =
    Then(Last)

Nth(N) = Query(Nth, N)

function First(env::Environment, p::Pipeline, X)
    x = assemble(env, target_pipe(p), X)
    compose(p, assemble_first(x))
end

function Last(env::Environment, p::Pipeline, X)
    x = assemble(env, target_pipe(p), X)
    compose(p, assemble_first(x, true))
end

function Nth(env::Environment, p::Pipeline, X, N::Int)
    x = assemble(env, target_pipe(p), X)
    compose(p, assemble_nth(x, N))
end

function Nth(env::Environment, p::Pipeline, X, N)
    p0 = target_pipe(p)
    x = assemble(env, p0, X)
    n = assemble(env, p0, N)
    compose(p, assemble_nth(x, n))
end

Nth(env::Environment, p::Pipeline, N::Int) =
    assemble_nth(p, N)

Nth(env::Environment, p::Pipeline, N) =
    assemble_nth(p, assemble(env, source_pipe(p), N))

translate(mod::Module, ::Val{:first}, (arg,)::Tuple{Any}) =
    First(translate(mod, arg))

translate(mod::Module, ::Val{:first}, ::Tuple{}) =
    Then(First)

translate(mod::Module, ::Val{:last}, (arg,)::Tuple{Any}) =
    Last(translate(mod, arg))

translate(mod::Module, ::Val{:last}, ::Tuple{}) =
    Then(Last)

translate(mod::Module, ::Val{:nth}, args::Tuple{Any,Any}) =
    Nth(translate.(Ref(mod), args)...)

translate(mod::Module, ::Val{:nth}, (arg,)::Tuple{Any}) =
    Nth(translate(mod, arg))

#
# Take and Drop combinators.
#

function assemble_take(p::Pipeline, n::Union{Int,Missing}, inv::Bool)
    elts = elements(target(p))
    card = cardinality(target(p))|x0to1
    chain_of(
        p,
        slice_by(n, inv),
    ) |> designate(source(p), BlockOf(elts, card) |> IsFlow)
end

function assemble_take(p::Pipeline, n::Pipeline, inv::Bool)
    n = uncover(n)
    fits(target(n), BlockOf(ValueOf(Int), x0to1)) || error("expected a singular integer")
    src = source(p)
    tgt = BlockOf(elements(target(p)), cardinality(target(p))|x0to1) |> IsFlow
    chain_of(
        tuple_of(p, n),
        slice_by(inv),
    ) |> designate(src, tgt)
end

"""
    Take(N) :: Query

This query preserves the first `N` elements of its input, dropping
the rest.

```jldoctest
julia> unitknot[Lift('a':'c') >> Take(2)]
──┼───┼
1 │ a │
2 │ b │
```

`Take(-N)` drops the last `N` elements.

```jldoctest
julia> unitknot[Lift('a':'c') >> Take(-2)]
──┼───┼
1 │ a │
```
"""
Take(N) =
    Query(Take, N)

"""
    Drop(N) :: Query

This query drops the first `N` elements of its input, preserving
the rest.

```jldoctest
julia> unitknot[Lift('a':'c') >> Drop(2)]
──┼───┼
1 │ c │
```

`Drop(-N)` takes the last `N` elements.

```jldoctest
julia> unitknot[Lift('a':'c') >> Drop(-2)]
──┼───┼
1 │ b │
2 │ c │
```
"""
Drop(N) =
    Query(Drop, N)

Take(env::Environment, p::Pipeline, n::Union{Int,Missing}, inv::Bool=false) =
    assemble_take(p, n, inv)

function Take(env::Environment, p::Pipeline, N, inv::Bool=false)
    n = assemble(env, source_pipe(p), N)
    assemble_take(p, n, inv)
end

Drop(env::Environment, p::Pipeline, N) =
    Take(env, p, N, true)

translate(mod::Module, ::Val{:take}, (arg,)::Tuple{Any}) =
    Take(translate(mod, arg))

translate(mod::Module, ::Val{:drop}, (arg,)::Tuple{Any}) =
    Drop(translate(mod, arg))

#
# Unique and Group combinators.
#

function assemble_unique(p::Pipeline, x::Pipeline)
    x = uncover(x)
    q = chain_of(x, unique_by()) |> designate(source(x), target(x)) |> cover
    compose(p, q)
end

"""
    Unique(X) :: Query

This query produces all distinct elements emitted by `X`.

```jldoctest
julia> unitknot[Unique(['a','b','b','c','c','c'])]
──┼───┼
1 │ a │
2 │ b │
3 │ c │
```

---

    Each(X >> Unique) :: Query

In the query form, `Unique` produces all distinct elements of its input.

```jldoctest
julia> unitknot[Lift(['a','b','b','c','c','c']) >> Unique]
──┼───┼
1 │ a │
2 │ b │
3 │ c │
```
"""
Unique(X) =
    Query(Unique, X)

Lift(::typeof(Unique)) =
    Then(Unique)

function Unique(env::Environment, p::Pipeline, X)
    x = assemble(env, target_pipe(p), X)
    assemble_unique(p, x)
end

translate(mod::Module, ::Val{:unique}, (arg,)::Tuple{Any}) =
    Unique(translate(mod, arg))

translate(mod::Module, ::Val{:unique}, ::Tuple{}) =
    Then(Unique)

function assemble_group(p::Pipeline, xs::Vector{Pipeline})
    lbls = Symbol[]
    ks = Pipeline[]
    seen = Dict{Symbol,Int}()
    for (i, x) in enumerate(xs)
        x = uncover(x)
        lbl = getlabel(x, nothing)
        if lbl !== nothing
            x = relabel(x, nothing)
        else
            lbl = ordinal_label(i)
        end
        if lbl in keys(seen)
            lbls[seen[lbl]] = ordinal_label(seen[lbl])
        end
        seen[lbl] = i
        push!(lbls, lbl)
        push!(ks, x)
    end
    p0 = uncover(target_pipe(p))
    lbl = getlabel(p0, nothing)
    if lbl !== nothing
        p0 = relabel(p0, nothing)
    else
        lbl = ordinal_label(length(xs)+1)
    end
    if lbl in keys(seen)
        lbls[seen[lbl]] = ordinal_label(seen[lbl])
    end
    push!(lbls, lbl)
    src = elements(target(p))
    tgt = TupleOf(target(p0), TupleOf(target.(ks)))
    q = tuple_of(p0, tuple_of(Symbol[], ks)) |> designate(src, tgt)
    r = uncover(compose(p, cover(q)))
    cols = Pipeline[]
    for (i, k) in enumerate(ks)
        col = chain_of(column(2), column(i)) |> designate(elements(target(r)), target(k))
        push!(cols, col)
    end
    card = cardinality(target(r)) & x1toN
    col = chain_of(column(1), flatten()) |> designate(elements(target(r)),
                                                      BlockOf(elements(target(p0)), card))
    push!(cols, col)
    src = source(r)
    tgt = replace_elements(target(r), TupleOf(lbls, target.(cols)))
    chain_of(
        r,
        group_by(),
        with_elements(tuple_of(lbls, cols)),
    ) |> designate(src, tgt) |> cover
end

"""
    Group(X₁, X₂ … Xₙ) :: Query

This query groups the input data by the keys `X₁`, `X₂` … `Xₙ`.

```jldoctest
julia> unitknot[Lift(1:5) >> Group(isodd.(It))]
  │ #A     #B      │
──┼────────────────┼
1 │ false  2; 4    │
2 │  true  1; 3; 5 │
```
"""
Group(Xs...) =
    Query(Group, Xs...)

function Group(env::Environment, p::Pipeline, Xs...)
    xs = assemble.(Ref(env), Ref(target_pipe(p)), collect(AbstractQuery, Xs))
    assemble_group(p, xs)
end

translate(mod::Module, ::Val{:group}, args::Tuple) =
    Group(translate.(Ref(mod), args)...)

#
# Cardinality assertions.
#

function assemble_cardinality(p::Pipeline, card::Cardinality)
    src = source(p)
    tgt = BlockOf(elements(target(p)), card) |> IsFlow
    q = block_cardinality(card, getlabel(src, nothing), getlabel(tgt, nothing))
    chain_of(p, q) |> designate(src, tgt)
end

"""
    Is0to1(X) :: Query

This query asserts that `X` emits 0 or 1 element.

---

    Each(X >> Is0to1) :: Query

In this form, `Is0to1` asserts that its input contains 0 or 1 element.
"""
Is0to1(X) = Query(Is0to1, X)

"""
    Is0toN(X) :: Query

This query asserts that `X` may emit any number of elements.

---

    Each(X >> Is0toN) :: Query

In this form, `Is0toN` asserts that its input contains any number of elements.
"""
Is0toN(X) = Query(Is0toN, X)

"""
    Is1to1(X) :: Query

This query asserts that `X` emits 1 element.

---

    Each(X >> Is1to1) :: Query

In this form, `Is1to1` asserts that its input contains 1 element.
"""
Is1to1(X) = Query(Is1to1, X)

"""
    Is1toN(X) :: Query

This query asserts that `X` emits 1 or more elements.

---

    Each(X >> Is1toN) :: Query

In this form, `Is1toN` asserts that its input contains 1 or more elements.
"""
Is1toN(X) = Query(Is1toN, X)

Lift(::typeof(Is0to1)) = Then(Is0to1)

Lift(::typeof(Is0toN)) = Then(Is0toN)

Lift(::typeof(Is1to1)) = Then(Is1to1)

Lift(::typeof(Is1toN)) = Then(Is1toN)

Is0to1(env::Environment, p::Pipeline, X) =
    assemble_cardinality(assemble(env, p, X), x0to1)

Is0toN(env::Environment, p::Pipeline, X) =
    assemble_cardinality(assemble(env, p, X), x0toN)

Is1to1(env::Environment, p::Pipeline, X) =
    assemble_cardinality(assemble(env, p, X), x1to1)

Is1toN(env::Environment, p::Pipeline, X) =
    assemble_cardinality(assemble(env, p, X), x1toN)

translate(mod::Module, ::Val{:is0to1}, (arg,)::Tuple{Any}) =
    Is0to1(translate(mod, arg))

translate(mod::Module, ::Val{:is0toN}, (arg,)::Tuple{Any}) =
    Is0toN(translate(mod, arg))

translate(mod::Module, ::Val{:is1to1}, (arg,)::Tuple{Any}) =
    Is1to1(translate(mod, arg))

translate(mod::Module, ::Val{:is1toN}, (arg,)::Tuple{Any}) =
    Is1toN(translate(mod, arg))

translate(mod::Module, ::Val{:is0to1}, ::Tuple{}) =
    Then(Is0to1)

translate(mod::Module, ::Val{:is0toN}, ::Tuple{}) =
    Then(Is0toN)

translate(mod::Module, ::Val{:is1to1}, ::Tuple{}) =
    Then(Is1to1)

translate(mod::Module, ::Val{:is1toN}, ::Tuple{}) =
    Then(Is1toN)


text:
# All Abstract types defined
"""
    AbstractKMeansAlg

Abstract base type inherited by all sub-KMeans algorithms.
"""
abstract type AbstractKMeansAlg end

"""
    ClusteringResult

Base type for the output of clustering algorithm.
"""
abstract type ClusteringResult end

# Here we mimic `Clustering` output structure
"""
    KmeansResult{C,D<:Real,WC<:Real} <: ClusteringResult

The output of [`kmeans`](@ref) and [`kmeans!`](@ref).
# Type parameters
 * `C<:AbstractMatrix{<:AbstractFloat}`: type of the `centers` matrix
 * `D<:Real`: type of the assignment cost
 * `WC<:Real`: type of the cluster weight
 # C is the type of centers, an (abstract) matrix of size (d x k)
# D is the type of pairwise distance computation from points to cluster centers
# WC is the type of cluster weights, either Int (in the case where points are
# unweighted) or eltype(weights) (in the case where points are weighted).
"""
struct KmeansResult{C<:AbstractMatrix{<:AbstractFloat},D<:Real,WC<:Real} <: ClusteringResult
    centers::C                 # cluster centers (d x k)
    assignments::Vector{Int}   # assignments (n)
    costs::Vector{D}           # cost of the assignments (n)
    counts::Vector{Int}        # number of points assigned to each cluster (k)
    wcounts::Vector{WC}        # cluster weights (k)
    totalcost::D               # total cost (i.e. objective)
    iterations::Int            # number of elapsed iterations
    converged::Bool            # whether the procedure converged
end

"""
    @parallelize(n_threads, ncol, f)

Parallelize function and run it over n_threads. Function should require following conditions:
1. It should not return any values.
1. It should accept parameters two parameters at the end of the argument list. First
accepted parameter is `range`, which defines chunk used in calculations. Second
parameter is `idx` which defines id of the container where results can be stored.

`ncol` argument defines range 1:ncol which is sliced in `n_threads` chunks.
"""
macro parallelize(n_threads, ncol, f)
    for i in 1:length(f.args)
        f.args[i] = :($(esc(f.args[i])))
    end
    single_thread_chunk = copy(f)
    push!(single_thread_chunk.args, :(1:$(esc(ncol))))
    push!(single_thread_chunk.args, 1)

    multi_thread_chunk = copy(f)
    push!(multi_thread_chunk.args, :(ranges[i]))
    push!(multi_thread_chunk.args, :(i))

    last_multi_thread_chunk = copy(f)
    push!(last_multi_thread_chunk.args, :(ranges[end]))
    push!(last_multi_thread_chunk.args, :($(esc(n_threads))))

    return quote
        if $(esc(n_threads)) == 1
            $single_thread_chunk
        else
            local ranges = splitter($(esc(ncol)), $(esc(n_threads)))
            local waiting_list = $(esc(Vector)){$(esc(Task))}(undef, $(esc(n_threads)) - 1)
            for i in 1:$(esc(n_threads)) - 1
                waiting_list[i] = @spawn $multi_thread_chunk
            end

            $last_multi_thread_chunk

            for i in 1:$(esc(n_threads)) - 1
                wait(waiting_list[i])
            end
        end
    end
end

"""
    distance(X1, X2, i1, i2)

Allocationless calculation of square eucledean distance between vectors X1[:, i1] and X2[:, i2]
"""
function distance(X1, X2, i1, i2)
    d = 0.0
    # TODO: break of the loop if d is larger than threshold (known minimum disatnce)
    @inbounds @simd for i in axes(X1, 1)
        d += (X1[i, i1] - X2[i, i2])^2
    end

    return d
end

"""
    sum_of_squares(x, labels, centre, k)

This function computes the total sum of squares based on the assigned (labels)
design matrix(x), centroids (centre), and the number of desired groups (k).

A Float type representing the computed metric is returned.
"""
function sum_of_squares(containers, x, labels, centre, r, idx)
    s = 0.0

    @inbounds for j in r
        for i in axes(x, 1)
            s += (x[i, j] - centre[i, labels[j]])^2
        end
    end

    containers.sum_of_squares[idx] = s
end

"""
    Kmeans([alg::AbstractKMeansAlg,] design_matrix, k; n_threads = nthreads(), k_init="k-means++", max_iters=300, tol=1e-6, verbose=true)

This main function employs the K-means algorithm to cluster all examples
in the training data (design_matrix) into k groups using either the
`k-means++` or random initialisation technique for selecting the initial
centroids.

At the end of the number of iterations specified (max_iters), convergence is
achieved if difference between the current and last cost objective is
less than the tolerance level (tol). An error is thrown if convergence fails.

Arguments:
- `alg` defines one of the algorithms used to calculate `k-means`. This
argument can be omitted, by default Lloyd algorithm is used.
- `n_threads` defines number of threads used for calculations, by default it is equal
to the `Threads.nthreads()` which is defined by `JULIA_NUM_THREADS` environmental
variable. For small size design matrices it make sense to set this argument to 1 in order
to avoid overhead of threads generation.
- `k_init` is one of the algorithms used for initialization. By default `k-means++` algorithm is used,
alternatively one can use `rand` to choose random points for init.
- `max_iters` is the maximum number of iterations
- `tol` defines tolerance for early stopping.
- `verbose` is verbosity level. Details of operations can be either printed or not by setting verbose accordingly.

A `KmeansResult` structure representing labels, centroids, and sum_squares is returned.
"""
function kmeans(alg, design_matrix, k;
                n_threads = Threads.nthreads(),
                k_init = "k-means++", max_iters = 300,
                tol = 1e-6, verbose = false, init = nothing)
    nrow, ncol = size(design_matrix)
    containers = create_containers(alg, k, nrow, ncol, n_threads)

    return kmeans!(alg, containers, design_matrix, k, n_threads = n_threads,
                    k_init = k_init, max_iters = max_iters, tol = tol,
                    verbose = verbose, init = init)
end


text:
# this is a bit of an odd design, but I think it does it's job.
# This file contains statements spereated by the comment # cell (without space),
# which makes it possible to extract each cell out of this file.
# There is also the setup cell indicated by # setup.
# We can then continue to generate files for visual regression tests
# snoop compiling and even create ijulia notebooks with cells containing this.

#setup
using Makie, GLFW, GeometryTypes, Reactive, FileIO, ColorBrewer, Colors
using GLVisualize
using GLVisualize: loadasset, assetpath

function xy_data(x, y)
    r = sqrt(x*x + y*y)
    r == 0.0 ? 1f0 : (sin(r)/r)
end

function custom_theme(scene)
    @theme theme = begin
        linewidth = to_float(3)
        colormap = to_colormap(:RdYlGn)#to_colormap(:RdPu)
        scatter = begin
            marker = to_spritemarker(Circle)
            markersize = to_float(0.03)
            strokecolor = to_color(:white)
            strokewidth = to_float(0.01)
            glowcolor = to_color(RGBA(0, 0, 0, 0.4))
            glowwidth = to_float(0.1)
        end
    end
    # update theme values
    scene[:theme] = theme
end

#cell
img = loadasset("doge.png")
scene = Scene(resolution = (500, 500))

display(scene)
show(scene)
println(scene)

is = image(img)
center!(scene)
subscene = Scene(scene, Signal(SimpleRectangle(0, 0, 200, 200)))
scatter(subscene, rand(100) * 200, rand(100) * 200, markersize = 4)
center!(scene)

#cell
scene = Scene(resolution = (500, 500));
x = [0, 1, 2, 0];
y = [0, 0, 1, 2];
z = [0, 2, 0, 1];
color = [:red, :green, :blue, :yellow];
i = [0, 0, 0, 1];
j = [1, 2, 3, 2];
k = [2, 3, 1, 3];

indices = [1, 2, 3, 1, 3, 4, 1, 4, 2, 2, 3, 4];
mesh(x, y, z, indices, color = color);
r = linspace(-0.5, 2.5, 4);
axis(r, r, r);
center!(scene);

#cell
scene = Scene(resolution = (500, 500))
Makie.Makie.volume(rand(32, 32, 32), algorithm = :iso)
center!(scene)

#cell
scene = Scene(resolution = (500, 500))
heatmap(rand(32, 32))
center!(scene)

#cell
scene = Scene(resolution = (500, 500))
r = linspace(-10, 10, 512)
z = ((x, y)-> sin(x) + cos(y)).(r, r')
Makie.contour(r, r, z, levels = 5, color = ColorBrewer.palette("RdYlBu", 5))
center!(scene)

#cell
scene = Scene(resolution = (500, 500))
vx = -1:0.1:1;
vy = -1:0.1:1;

f(x, y) = (sin(x*10) + cos(y*10)) / 4
psurf = surface(vx, vy, f)

pos = lift_node(psurf[:x], psurf[:y], psurf[:z]) do x, y, z
    vec(Point3f0.(x, y', z .+ 0.5))
end
pscat = scatter(pos)
plines = lines(view(pos, 1:2:length(pos)))
center!(scene)
@theme theme = begin
    markersize = to_markersize2d(0.01)
    strokecolor = to_color(:white)
    strokewidth = to_float(0.01)
end
# this pushes all the values from theme to the plot
push!(pscat, theme)
pscat[:glow_color] = to_node(RGBA(0, 0, 0, 0.4), x->to_color((), x))
# apply it to the scene
custom_theme(scene)
# From now everything will be plotted with new theme
psurf = surface(vx, 1:0.1:2, psurf[:z])
center!(scene)

#cell
scene = Scene(resolution = (500, 500))
sv = scatter(rand(Point3f0, 100))
similar(sv, rand(10), rand(10), rand(10), color = :black, markersize = 0.4)
scene

#cell
scene = Scene(resolution = (500, 500))
large_sphere = HyperSphere(Point3f0(0), 1f0)
positions = decompose(Point3f0, large_sphere)
colS = [Colors.RGBA{Float32}(rand(), rand(), rand(), 1.) for i = 1:length(positions)]
sizesS = [rand(Vec3f0) .* 0.5f0 for i = 1:length(positions)]
meshscatter(positions, color = colS, markersize = sizesS)
center!(scene)

#cell
scene = Scene(resolution = (500, 500))
y = [
    -0.997669
    -0.979084
    -0.942261
    -0.887885
    -0.81697
    -0.730836
    -0.631088
    -0.519584
    -0.398401
    -0.269797
    -0.136167
    0.0
    0.136167
    0.269797
    0.398401
    0.519584
    0.631088
    0.730836
    0.81697
    0.887885
    0.942261
    0.979084
    0.997669
]
contour(linspace(-0.99, 0.99, 23), y, rand(23, 23), levels = 10)
center!(scene)

#cell
scene = Scene(resolution = (500, 500))
# define points/edges
perturbfactor = 4e1
N = 3; nbfacese = 30; radius = 0.02
large_sphere = HyperSphere(Point3f0(0), 1f0)
positions = decompose(Point3f0, large_sphere, 30)
np = length(positions)
pts = [positions[k][l] for k = 1:length(positions), l = 1:3]
pts = vcat(pts, 1.1 * pts + randn(size(pts)) / perturbfactor) # light position influence ?
edges = hcat(collect(1:np), collect(1:np) + np)
ne = size(edges, 1); np = size(pts, 1)
# define markers meshes
meshC = GeometryTypes.GLNormalMesh(GeometryTypes.Cylinder{3, Float32}(
                                   GeometryTypes.Point3f0(0., 0., 0.),
                                   GeometryTypes.Point3f0(0., 0, 1.),
                                   Float32(1)), nbfacese)

meshS = GeometryTypes.GLNormalMesh(large_sphere, 20)
# define colors, markersizes and rotations
pG = [GeometryTypes.Point3f0(pts[k, 1], pts[k, 2], pts[k, 3]) for k = 1:np]
lengthsC = sqrt.(sum((pts[edges[:,1], :] .- pts[edges[:, 2], :]) .^ 2, 2))
sizesC = [GeometryTypes.Vec3f0(radius, radius, lengthsC[i]) for i = 1:ne]
sizesC = [Vec3f0(1., 1., 1.) for i = 1:ne]
colorsp = [Colors.RGBA{Float32}(rand(), rand(), rand(), 1.) for i = 1:np]
colorsC = [(colorsp[edges[i, 1]] + colorsp[edges[i, 2]]) / 2. for i = 1:ne]
sizesC = [Vec3f0(radius, radius, lengthsC[i]) for i = 1:ne]
Qlist = zeros(ne, 4)
for k = 1:ne
    ct = GeometryTypes.Cylinder{3, Float32}(
                GeometryTypes.Point3f0(pts[edges[k, 1], 1], pts[edges[k, 1], 2], pts[edges[k, 1], 3]),
                GeometryTypes.Point3f0(pts[edges[k, 2], 1], pts[edges[k, 2], 2], pts[edges[k, 2], 3]),
                Float32(1))
    Q = GeometryTypes.rotation(ct)
    r = 0.5 * sqrt(1 + Q[1, 1] + Q[2, 2] + Q[3, 3]); Qlist[k, 4] = r
    Qlist[k, 1] = (Q[3, 2] - Q[2, 3]) / (4 * r)
    Qlist[k, 2] = (Q[1, 3] - Q[3, 1]) / (4 * r)
    Qlist[k, 3] = (Q[2, 1] - Q[1, 2]) / (4 * r)
end
rotationsC = AbstractVector[Vec4f0(Qlist[i, 1], Qlist[i, 2], Qlist[i, 3], Qlist[i, 4]) for i = 1:ne]
# plot
hm = Makie.meshscatter(pG[edges[:, 1]], color = colorsC, marker = meshC,
                       markersize = sizesC,  rotations = rotationsC)
hp = Makie.meshscatter(pG, color = colorsp, marker = meshS, markersize = radius)

r = linspace(-1.3, 1.3, 4); Makie.axis(r, r, r)
center!(scene)

#cell
scene = Scene(resolution = (500, 500))
large_sphere = HyperSphere(Point3f0(0), 1f0)
positions = decompose(Point3f0, large_sphere)
linepos = view(positions, rand(1:length(positions), 1000))
lines(linepos, linewidth = 0.1, color = :black)
scatter(positions, strokewidth = 0.02, strokecolor = :white, color = RGBA(0.9, 0.2, 0.4, 0.6))
r = linspace(-1.5, 1.5, 5)
axis(r, r, r)
scene

#cell
scene = Scene(resolution = (500, 500))
large_sphere = HyperSphere(Point3f0(0), 1f0)
positions = decompose(Point3f0, large_sphere)
meshscatter(positions, color = RGBA(0.9, 0.2, 0.4, 1))
scene

#cell
scene = Scene(resolution = (500, 500))

r = linspace(-2, 2, 40)
surf_func(i) = [Float32(xy_data(x*i, y*i)) for x = r, y = r]
z = surf_func(20)
surf = surface(r, r, z)

wf = wireframe(r, r, surf[:z] .+ 1.0,
    linewidth = 2f0, color = lift_node(x-> x[5], surf[:colormap])
)
xy = linspace(-2.1, 2.1, 4)
axis(xy, xy, linspace(0, 2, 4))
center!(scene)

io = VideoStream(scene)
for i in linspace(0, 60, 100)
    surf[:z] = surf_func(i)
    recordframe!(io)
end
scene

#cell
scene = Scene(resolution = (500, 500))

N = 40
r = linspace(-2, 2, 40)
surf_func(i) = [Float32(xy_data(x*i, y*i)) for x = r, y = r]
surface(
    r, r, surf_func(10),
    color = GLVisualize.loadasset("doge.png")
)
center!(scene)
scene

#cell
scene = Scene(resolution = (500, 500))
x = GLVisualize.loadasset("cat.obj")
Makie.mesh(x.vertices, x.faces, color = :black)
pos = map(x.vertices, x.normals) do p, n
    p => p .+ (normalize(n) .* 0.05f0)
end
linesegment(pos)
scene

#cell
scene = Scene(resolution = (500, 500))
mesh(GLVisualize.loadasset("cat.obj"))
r = linspace(-0.1, 1, 4)
center!(scene)
scene

#cell
scene = Scene(resolution = (500, 500))
cat = load(assetpath("cat.obj"), GLNormalUVMesh)
Makie.mesh(cat, color = loadasset("diffusemap.tga"))
center!(scene)

#cell
scene = Scene(resolution = (500, 500))
Makie.mesh(Sphere(Point3f0(0), 1f0))
center!(scene)
scene

#cell
scene = Scene(resolution = (500, 500))
wireframe(GLVisualize.loadasset("cat.obj"))
center!(scene)
scene

#cell
scene = Scene(resolution = (500, 500))
wireframe(Sphere(Point3f0(0), 1f0))
center!(scene)
scene

#cell
scene = Scene(resolution = (500, 500))
heatmap(rand(32, 32))
center!(scene)

#cell
scene = Scene(resolution = (500, 500), color = :black)
stars = 100_000
scatter((rand(Point3f0, stars) .- 0.5) .* 10,
    glowwidth = 0.005, glow_color = :white, color = RGBA(0.8, 0.9, 0.95, 0.4),
    markersize = rand(linspace(0.0001, 0.01, 100), stars)
)
scene

#cell
scene = Scene(resolution = (500, 500))
Makie.volume(rand(32, 32, 32), algorithm = :iso)
center!(scene)

#cell
scene = Scene(resolution = (500, 500))
scatter(Point3f0[(1,0,0), (0,1,0), (0,0,1)], marker=[:x, :circle, :cross])
axis(scene, linspace(0, 1, 4), linspace(0, 1, 4), linspace(0, 1, 4))
center!(scene);

#cell
scene = Scene(resolution = (500, 500))

x = map([:dot, :dash, :dashdot], [2, 3, 4]) do ls, lw
    linesegment(linspace(1, 5, 100), rand(100), rand(100), linestyle = ls, linewidth = lw)
end
push!(x, scatter(linspace(1, 5, 100), rand(100), rand(100)))
center!(scene)
l = Makie.legend(x, ["attribute $i" for i in 1:4])
l[:position] = (0, 1)
l[:backgroundcolor] = RGBA(0.95, 0.95, 0.95)
l[:strokecolor] = RGB(0.8, 0.8, 0.8)
l[:gap] = 30
l[:textsize] = 19
l[:linepattern] = Point2f0[(0,-0.2), (0.5, 0.2), (0.5, 0.2), (1.0, -0.2)]
l[:scatterpattern] = decompose(Point2f0, Circle(Point2f0(0.5, 0), 0.3f0), 9)
l[:markersize] = 2f0
scene

#cell
scene = Scene(resolution = (500, 500))
cmap = collect(linspace(to_color(:red), to_color(:blue), 20))
l = Makie.legend(cmap, 1:4)
l[:position] = (1.0,1.0)
l[:textcolor] = :blue
l[:strokecolor] = :black
l[:strokewidth] = 1
l[:textsize] = 15
l[:textgap] = 5
scene