text:
abstract type AbstractBSONReader end

struct BSONReader{S <: DenseVector{UInt8}, V <: BSONValidator} <: AbstractBSONReader
    src::S
    offset::Int
    type::UInt8
    validator::V
end

@inline function BSONReader(src::DenseVector{UInt8}, validator = LightBSONValidator())
    validate_root(validator, src)
    BSONReader(src, 0, BSON_TYPE_DOCUMENT, validator)
end

@inline Base.pointer(reader::BSONReader) = pointer(reader.src) + reader.offset

const TYPE_SIZE_TABLE = fill(-1, 256)
TYPE_SIZE_TABLE[BSON_TYPE_DOUBLE] = 8
TYPE_SIZE_TABLE[BSON_TYPE_DATETIME] = 8
TYPE_SIZE_TABLE[BSON_TYPE_INT64] = 8
TYPE_SIZE_TABLE[BSON_TYPE_TIMESTAMP] = 8
TYPE_SIZE_TABLE[BSON_TYPE_INT32] = 4
TYPE_SIZE_TABLE[BSON_TYPE_BOOL] = 1
TYPE_SIZE_TABLE[BSON_TYPE_NULL] = 0
TYPE_SIZE_TABLE[BSON_TYPE_UNDEFINED] = 0
TYPE_SIZE_TABLE[BSON_TYPE_DECIMAL128] = 16
TYPE_SIZE_TABLE[BSON_TYPE_OBJECTID] = 12
TYPE_SIZE_TABLE[BSON_TYPE_MIN_KEY] = 0
TYPE_SIZE_TABLE[BSON_TYPE_MAX_KEY] = 0

function element_size_variable_(t::UInt8, p::Ptr{UInt8})
    if t == BSON_TYPE_DOCUMENT || t == BSON_TYPE_ARRAY || t == BSON_TYPE_CODE_WITH_SCOPE
        Int(ltoh(unsafe_load(Ptr{Int32}(p))))
    elseif t == BSON_TYPE_STRING || t == BSON_TYPE_CODE || t == BSON_TYPE_SYMBOL
        Int(ltoh(unsafe_load(Ptr{Int32}(p)))) + 4
    elseif t == BSON_TYPE_BINARY
        Int(ltoh(unsafe_load(Ptr{Int32}(p)))) + 5
    elseif t == BSON_TYPE_REGEX
        len1 = unsafe_trunc(Int, ccall(:strlen, Csize_t, (Cstring,), p)) + 1
        len2 = unsafe_trunc(Int, ccall(:strlen, Csize_t, (Cstring,), p + len1)) + 1
        len1 + len2
    elseif t == BSON_TYPE_DB_POINTER
        Int(ltoh(unsafe_load(Ptr{Int32}(p)))) + 16
    else
        error("Unsupported BSON type $t")
    end
end

@inline function element_size_(t::UInt8, p::Ptr{UInt8})
    @inbounds s = TYPE_SIZE_TABLE[t]
    s >= 0 ? s : element_size_variable_(t, p)
end

"""
    name_len_and_match_(field, target)

Compare field name and target name, and calculate length of field name (including null terminator)

Returns (length(field)), field == target)

Performs byte wise comparison, optimized for short length field names
"""
@inline function name_len_and_match_(field::Ptr{UInt8}, target::Ptr{UInt8})
    len = 1
    cmp_acc = 0x0
    while true
        c = unsafe_load(field, len)
        c2 = unsafe_load(target, len)
        cmp_acc |= xor(c, c2)
        c == 0x0 && return (len, cmp_acc == 0)
        len += 1
        c2 == 0x0 && break
    end
    while true
        c = unsafe_load(field, len)
        c == 0x0 && break
        len += 1
    end
    (len, false)
end

@inline function name_len_(field::Ptr{UInt8})
    len = 1
    while true
        unsafe_load(field, len) == 0x0 && break
        len += 1
    end
    len
end

@inline function Transducers.__foldl__(rf, val, reader::BSONReader)
    reader.type == BSON_TYPE_DOCUMENT || reader.type == BSON_TYPE_ARRAY || throw(
        ArgumentError("Field access only available on documents and arrays")
    )
    src = reader.src
    GC.@preserve src begin
        p = pointer(reader.src)
        offset = reader.offset
        doc_len = Int(ltoh(unsafe_load(Ptr{Int32}(p + offset))))
        doc_end = offset + doc_len - 1
        offset += 4
        while offset < doc_end
            el_type = unsafe_load(p, offset + 1)
            name_p = p + offset + 1
            name_len = name_len_(name_p)
            field_p = name_p + name_len
            value_len = element_size_(el_type, field_p)
            field_start = offset + 1 + name_len
            field_end = field_start + value_len
            validate_field(reader.validator, el_type, field_p, value_len, doc_end - field_start)
            field_reader = BSONReader(src, field_start, el_type, reader.validator)
            val = Transducers.@next(rf, val, UnsafeBSONString(name_p, name_len - 1) => field_reader)
            offset = field_end
        end
        Transducers.complete(rf, val)
    end
end

@inline function Base.foreach(f, reader::BSONReader)
    foreach(f, Map(identity), reader)
end

function Base.getindex(reader::BSONReader, target::Union{AbstractString, Symbol})
    reader.type == BSON_TYPE_DOCUMENT || reader.type == BSON_TYPE_ARRAY || throw(
        ArgumentError("Field access only available on documents and arrays")
    )
    src = reader.src
    GC.@preserve src target begin
        target_p = Base.unsafe_convert(Ptr{UInt8}, target)
        p = pointer(reader.src)
        offset = reader.offset
        doc_len = Int(ltoh(unsafe_load(Ptr{Int32}(p + offset))))
        doc_end = offset + doc_len - 1
        offset += 4
        while offset < doc_end
            el_type = unsafe_load(p, offset + 1)
            name_p = p + offset + 1
            name_len, name_match = name_len_and_match_(name_p, target_p)
            field_p = name_p + name_len
            value_len = element_size_(el_type, field_p)
            field_start = offset + 1 + name_len
            field_end = field_start + value_len
            validate_field(reader.validator, el_type, field_p, value_len, doc_end - field_start)
            name_match && return BSONReader(reader.src, field_start, el_type, reader.validator)
            offset = field_end
        end
    end
    throw(KeyError(target))
end

function Base.getindex(reader::BSONReader, i::Integer)
    i < 1 && throw(BoundsError(reader, i))
    el = foldl((_, x) -> reduced(x.second), Drop(i - 1), reader; init = nothing)
    el === nothing && throw(BoundsError(reader, i))
    el
end

@inline load_bits_(::Type{T}, p::Ptr{UInt8}) where T = ltoh(unsafe_load(Ptr{T}(p)))
@inline load_bits_(::Type{BSONTimestamp}, p::Ptr{UInt8}) = BSONTimestamp(load_bits_(UInt64, p))
@inline load_bits_(::Type{BSONObjectId}, p::Ptr{UInt8}) = unsafe_load(Ptr{BSONObjectId}(p))
@inline load_bits_(::Type{DateTime}, p::Ptr{UInt8}) = DateTime(
    Dates.UTM(Dates.UNIXEPOCH + load_bits_(Int64, p))
)

@inline function read_field_(reader::BSONReader, ::Type{T}) where T <: Union{
    Int64,
    Int32,
    Float64,
    Dec128,
    DateTime,
    BSONTimestamp,
    BSONObjectId
}
    reader.type == bson_type_(T) || throw(BSONConversionError(reader.type, T))
    GC.@preserve reader load_bits_(T, pointer(reader))
end

@inline function read_field_(reader::BSONReader, ::Type{Bool})
    reader.type == BSON_TYPE_BOOL || throw(BSONConversionError(reader.type, Bool))
    x = GC.@preserve reader unsafe_load(pointer(reader))
    validate_bool(reader.validator, x)
    x != 0x0
end

@inline function read_binary_(reader::BSONReader)
    GC.@preserve reader begin
        p = pointer(reader)
        len = load_bits_(Int32, p)
        subtype = unsafe_load(p + 4)
        validate_binary_subtype(reader.validator, p, len, subtype)
        UnsafeBSONBinary(
            UnsafeArray(p + 5, (Int(len),)),
            subtype
        )
    end
end

@inline function read_field_(reader::BSONReader, ::Type{UnsafeBSONBinary})
    reader.type == BSON_TYPE_BINARY || throw(BSONConversionError(reader.type, UnsafeBSONBinary))
    read_binary_(reader)
end

function read_field_(reader::BSONReader, ::Type{BSONBinary})
    reader.type == BSON_TYPE_BINARY || throw(BSONConversionError(reader.type, BSONBinary))
    x = read_binary_(reader)
    GC.@preserve reader BSONBinary(copy(x.data), x.subtype)
end

@inline function read_field_(reader::BSONReader, ::Type{UUID})
    reader.type == BSON_TYPE_BINARY || throw(BSONConversionError(reader.type, UUID))
    GC.@preserve reader begin
        p = pointer(reader)
        subtype = unsafe_load(p + 4)
        subtype == BSON_SUBTYPE_UUID || subtype == BSON_SUBTYPE_UUID_OLD || throw(
            BSONConversionError(reader.type, subtype, UUID)
        )
        len = load_bits_(Int32, p)
        len != 16 && error("Unexpected UUID length $len")
        return unsafe_load(Ptr{UUID}(p + 5))
    end
end

@inline function read_field_(reader::BSONReader, ::Type{BSONUUIDOld})
    reader.type == BSON_TYPE_BINARY || throw(BSONConversionError(reader.type, BSONUUIDOld))
    GC.@preserve reader begin
        p = pointer(reader)
        subtype = unsafe_load(p + 4)
        subtype == BSON_SUBTYPE_UUID_OLD || throw(BSONConversionError(reader.type, subtype, BSONUUIDOld))
        len = load_bits_(Int32, p)
        len != 16 && error("Unexpected UUID length $len")
        return BSONUUIDOld(unsafe_load(Ptr{UUID}(p + 5)))
    end
end

@inline function read_string_(reader::BSONReader)
    GC.@preserve reader begin
        p = pointer(reader)
        len = Int(load_bits_(Int32, p))
        validate_string(reader.validator, p + 4, len - 1)
        UnsafeBSONString(p + 4, len - 1)
    end
end

@inline function read_field_(reader::BSONReader, ::Type{UnsafeBSONString})
    reader.type == BSON_TYPE_STRING || reader.type == BSON_TYPE_CODE || reader.type == BSON_TYPE_SYMBOL || throw(
        BSONConversionError(reader.type, UnsafeBSONString)
    )
    read_string_(reader)
end

function read_field_(reader::BSONReader, ::Type{String})
    reader.type == BSON_TYPE_STRING || reader.type == BSON_TYPE_CODE || reader.type == BSON_TYPE_SYMBOL || throw(
        BSONConversionError(reader.type, String)
    )
    GC.@preserve reader String(read_string_(reader))
end

function read_field_(reader::BSONReader, ::Type{BSONCode})
    reader.type == BSON_TYPE_CODE || throw(BSONConversionError(reader.type, BSONCode))
    GC.@preserve reader begin
        p = pointer(reader)
        len = Int(load_bits_(Int32, p))
        validate_string(reader.validator, p + 4, len - 1)
        BSONCode(unsafe_string(p + 4, len - 1))
    end
end

function read_field_(reader::BSONReader, ::Type{BSONSymbol})
    reader.type == BSON_TYPE_SYMBOL || throw(BSONConversionError(reader.type, BSONSymbol))
    BSONSymbol(reader[String])
end

function read_field_(reader::BSONReader, ::Type{BSONRegex})
    reader.type == BSON_TYPE_REGEX || throw(BSONConversionError(reader.type, BSONRegex))
    GC.@preserve reader begin
        p = pointer(reader)
        len1 = unsafe_trunc(Int, ccall(:strlen, Csize_t, (Cstring,), p))
        BSONRegex(
            unsafe_string(p, len1),
            unsafe_string(p + len1 + 1)
        )
    end
end

function read_field_(reader::BSONReader, ::Type{BSONDBPointer})
    reader.type == BSON_TYPE_DB_POINTER || throw(BSONConversionError(reader.type, BSONDBPointer))
    GC.@preserve reader begin
        p = pointer(reader)
        len = Int(load_bits_(Int32, p))
        validate_string(reader.validator, p + 4, len - 1)
        collection = unsafe_string(p + 4, len - 1)
        ref = unsafe_load(Ptr{BSONObjectId}(p + 4 + len))
        BSONDBPointer(collection, ref)
    end
end

@inline function read_field_(reader::BSONReader, ::Type{BSONMinKey})
    reader.type == BSON_TYPE_MIN_KEY || throw(BSONConversionError(reader.type, BSONMinKey))
    BSONMinKey()
end

@inline function read_field_(reader::BSONReader, ::Type{BSONMaxKey})
    reader.type == BSON_TYPE_MAX_KEY || throw(BSONConversionError(reader.type, BSONMaxKey))
    BSONMaxKey()
end

@inline function read_field_(reader::BSONReader, ::Type{BSONUndefined})
    reader.type == BSON_TYPE_UNDEFINED || throw(BSONConversionError(reader.type, BSONUndefined))
    BSONUndefined()
end

@inline function read_field_(reader::BSONReader, ::Type{Nothing})
    reader.type == BSON_TYPE_NULL || throw(BSONConversionError(reader.type, Nothing))
    nothing
end

function read_field_(reader::AbstractBSONReader, ::Type{Number})
    t = reader.type
    if t == BSON_TYPE_DOUBLE
        reader[Float64]
    elseif t == BSON_TYPE_INT64
        reader[Int64]
    elseif t == BSON_TYPE_INT32
        reader[Int32]
    elseif t == BSON_TYPE_DECIMAL128
        reader[Dec128]
    else
        throw(BSONConversionError(t, Number))
    end
end

function read_field_(reader::AbstractBSONReader, ::Type{Integer})
    t = reader.type
    if t == BSON_TYPE_INT64
        reader[Int64]
    elseif t == BSON_TYPE_INT32
        reader[Int32]
    else
        throw(BSONConversionError(t, Number))
    end
end

function read_field_(reader::AbstractBSONReader, ::Type{AbstractFloat})
    t = reader.type
    if t == BSON_TYPE_DOUBLE
        reader[Float64]
    elseif t == BSON_TYPE_DECIMAL128
        reader[Dec128]
    else
        throw(BSONConversionError(t, Number))
    end
end

@inline function read_field_(reader::AbstractBSONReader, ::Type{Union{Nothing, T}}) where T
    reader.type == BSON_TYPE_NULL ? nothing : reader[T]
end

function read_field_(reader::BSONReader, ::Type{BSONCodeWithScope})
    reader.type == BSON_TYPE_CODE_WITH_SCOPE || throw(BSONConversionError(reader.type, BSONCodeWithScope))
    src = reader.src
    GC.@preserve src begin
        offset = reader.offset
        p = pointer(reader.src) + offset
        code_len = Int(load_bits_(Int32, p + 4))
        doc_offset = offset + code_len + 8
        validate_field(reader.validator, BSON_TYPE_CODE, p + 8, code_len, Int(load_bits_(Int32, p)) - 11)
        validate_string(reader.validator, p + 8, code_len - 1)
        BSONCodeWithScope(
            unsafe_string(p + 8, code_len - 1),
            BSONReader(src, doc_offset, BSON_TYPE_DOCUMENT, reader.validator)[Dict{String, Any}]
        )
    end
end

function read_field_(reader::AbstractBSONReader, ::Type{T}) where T <: AbstractDict{String, Any}
    foldxl(reader; init = T()) do state, x
        state[String(x.first)] = x.second[Any]
        state
    end
end

function read_field_(reader::AbstractBSONReader, ::Type{T}) where T <: AbstractDict{Symbol, Any}
    foldxl(reader; init = T()) do state, x
        state[Symbol(x.first)] = x.second[Any]
        state
    end
end

function read_field_(reader::AbstractBSONReader, ::Type{Vector{T}}) where T
    dst = T[]
    copy!(dst, reader)
end

function Base.copy!(dst::AbstractArray{T}, reader::AbstractBSONReader) where T
    copy!(Map(x -> x.second[T]), dst, reader)
end

function read_field_(reader::AbstractBSONReader, ::Type{Any})
    if reader.type == BSON_TYPE_DOUBLE
        reader[Float64]
    elseif reader.type == BSON_TYPE_STRING
        reader[String]
    elseif reader.type == BSON_TYPE_DOCUMENT
        reader[OrderedDict{String, Any}]
    elseif reader.type == BSON_TYPE_ARRAY
        reader[Vector{Any}]
    elseif reader.type == BSON_TYPE_BINARY
        x = reader[BSONBinary]
        data = x.data
        if x.subtype == BSON_SUBTYPE_UUID && length(data) == 16
            GC.@preserve data unsafe_load(Ptr{UUID}(pointer(data)))
        elseif x.subtype == BSON_SUBTYPE_UUID_OLD && length(data) == 16
            GC.@preserve data BSONUUIDOld(unsafe_load(Ptr{UUID}(pointer(data))))
        else
            x
        end
    elseif reader.type == BSON_TYPE_OBJECTID
        reader[BSONObjectId]
    elseif reader.type == BSON_TYPE_BOOL
        reader[Bool]
    elseif reader.type == BSON_TYPE_DATETIME
        reader[DateTime]
    elseif reader.type == BSON_TYPE_NULL
        nothing
    elseif reader.type == BSON_TYPE_REGEX
        reader[BSONRegex]
    elseif reader.type == BSON_TYPE_CODE
        reader[BSONCode]
    elseif reader.type == BSON_TYPE_SYMBOL
        reader[BSONSymbol]
    elseif reader.type == BSON_TYPE_INT32
        reader[Int32]
    elseif reader.type == BSON_TYPE_TIMESTAMP
        reader[BSONTimestamp]
    elseif reader.type == BSON_TYPE_INT64
        reader[Int64]
    elseif reader.type == BSON_TYPE_DECIMAL128
        reader[Dec128]
    elseif reader.type == BSON_TYPE_MIN_KEY
        reader[BSONMinKey]
    elseif reader.type == BSON_TYPE_MAX_KEY
        reader[BSONMaxKey]
    elseif reader.type == BSON_TYPE_UNDEFINED
        reader[BSONUndefined]
    elseif reader.type == BSON_TYPE_CODE_WITH_SCOPE
        reader[BSONCodeWithScope]
    elseif reader.type == BSON_TYPE_DB_POINTER
        reader[BSONDBPointer]
    else
        error("Unsupported BSON type $(reader.type)")
    end
end

@inline @generated function bson_read_simple(::Type{T}, reader::AbstractBSONReader) where T
    field_readers = map(zip(fieldnames(T), fieldtypes(T))) do (fn, ft)
        fns = string(fn)
        :(reader[$fns][$ft])
    end
    :($T($(field_readers...)))
end

@inline function bson_read_structtype(::Type{T}, reader::AbstractBSONReader) where T
    StructTypes.construct((i, name, FT) -> reader[name][FT], T)
end

@inline function bson_read_unversioned(::Type{T}, reader::AbstractBSONReader) where T
    if bson_simple(T)
        bson_read_simple(T, reader)
    else
        bson_read_structtype(T, reader)
    end
end

@inline function bson_read_versioned(::Type{T}, v::V, reader::AbstractBSONReader) where {T, V}
    target = bson_schema_version(T)
    v != target && error("Mismatched schema version, read: $v, target: $target")
    bson_read_unversioned(T, reader)
end

@inline function bson_read(::Type{T}, reader::AbstractBSONReader) where T
    target_v = bson_schema_version(T)
    if target_v !== nothing
        v = reader[bson_schema_version_field(T)][typeof(target_v)]
        bson_read_versioned(T, v, reader)
    else
        bson_read_unversioned(T, reader)
    end
end

@inline function read_field_(reader::AbstractBSONReader, ::Type{T}) where T
    if isstructtype(T)
        bson_read(T, reader)
    else
        throw(ArgumentError("Unsupported type $T"))
    end
end

@inline function Base.getindex(reader::AbstractBSONReader, ::Type{T}) where T
    RT = bson_representation_type(T)
    if RT != T
        bson_representation_convert(T, read_field_(reader, RT))
    else
        read_field_(reader, T)
    end
end

@inline @generated function Base.getindex(reader::AbstractBSONReader, ::Type{T}) where T <: NamedTuple
    field_readers = map(zip(fieldnames(T), fieldtypes(T))) do (fn, ft)
        fns = string(fn)
        :(reader[$fns][$ft])
    end
    :($T(($(field_readers...),)))
end


text:
import Base.copy

export
    Graph,
    add_node!,
    add_undirected!

"""
    struct Graph

A simple graph with nodes (Int), edges (outgoing), and a size
property for each node
"""
struct Graph
    nodes :: Array{Int, 1}
    edges :: Dict{Int, Array{Int, 1}}
    sizes :: Dict{Int, Int}
    Graph() = new([], Dict(), Dict())
    Graph(ns, es, ss) = new(ns, es, ss)
end

function add_node!(g :: Graph, n :: Int, size :: Int)
    if !(n in g.nodes)
        push!(g.nodes, n)
        g.edges[n] = []
    end
    g.sizes[n] = size
end

function add_undirected!(g :: Graph, n1 :: Int, n2 :: Int)
    if n1 in g.nodes && n2 in g.nodes
        if !(n2 in g.edges[n1])
            push!(g.edges[n1], n2)
        end
        if !(n1 in g.edges[n2])
            push!(g.edges[n2], n1)
        end
    end
end

ext: jl


text:
### Utility Functions

"""
symarr_to_symengine(symarr::Vector{Symbol})

Converts a Vector{Symbol} into a string for SymEngine parsing

Symbol[:x,:y] --> "x,y"
"""
function symarr_to_symengine(symarr::Vector{Symbol})
  str = ""
  for sym in symarr
    str = str*string(sym)*","
  end
  str[1:end-1]
end


text:
function similar_normalized_eye(C::T) where {T<:AbstractMatrix}
    @assert size(C, 1) == size(C, 2)
    T(Diagonal(I / sqrt(size(C, 1)), size(C, 1)))
end

function _tensor_to_matrix(X::AbstractTensor3, lp_r::Bool = true)
    d = size(X, 1)
    lp_r ? reshape(permutedims(X, (1, 3, 2)), :, d) : reshape(X, d, :)
end

function _matrix_to_tensor(X::AbstractMatrix, lp_r::Bool = true)
    d = lp_r ? size(X, 2) : size(X, 1)
    lp_r ? permutedims(reshape(X, d, :, d), (1, 3, 2)) : reshape(X, d, d, :)
end

"""
    al_from_ac_and_c(AC, C)

Compute AL which satisfies `mul_matrix_from_right(AL, C) ≈ AC`

### Return values:
`AL, ϵL = al_from_ac_and_c(AC, C)`

where
`ϵL = norm(mul_matrix_from_right(AL, C) - AC)`
"""
function al_from_ac_and_c(AC::AbstractTensor3, C::AbstractMatrix)
    AC_Cdag = mul_matrix_from_right(AC, C')
    x, y, z = svd(_tensor_to_matrix(AC_Cdag, true))
    AL = _matrix_to_tensor(x * z', true)
    ϵL = norm(AC - mul_matrix_from_right(AL, C))

    AL, ϵL
end

"""
    ar_from_ac_and_c(AC, C)

Compute AR which satisfies `mul_matrix_from_left(AR, C) ≈ AC`

### Return values:
`AR, ϵR = ar_from_ac_and_c(AC, C)`

where
`ϵR = norm(mul_matrix_from_left(AR, C) - AC)`
"""
function ar_from_ac_and_c(AC::AbstractTensor3, C::AbstractMatrix)
    Cdag_AC = mul_matrix_from_left(AC, C')
    x, y, z = svd(_tensor_to_matrix(Cdag_AC, false))
    AR = _matrix_to_tensor(x * z', false)
    ϵR = norm(AC - mul_matrix_from_left(AR, C))

    AR, ϵR
end

"""
    mixedcanonical(A)

convert non-canonical uniform MPS `|Ψ(A)⟩ = ∑ vₗ^† (∏ᵢ A^{sᵢ}) vᵣ |{s}⟩` to
mixed-canonical uniform MPS
`|Ψ(AL, AR, AC, C)⟩ = ∑ vₗ^† (∏ᵢ AL^{sᵢ}) AC^{sⱼ} (∏ᵢ AR^{sᵢ}) vᵣ |{s}⟩ = ∑ vₗ^† (∏ᵢ AL^{sᵢ}) C (∏ᵢ AR^{sᵢ}) vᵣ |{s}⟩`

### Return values:
`AL, AR, AC, C = mixedcanonical(A)`
"""
function mixedcanonical(A::Tensor3{T}) where {T}

    d = size(A, 1)

    (L_val,), (L_vec,), L_inf = eigsolve(randn(T, d, d)) do X
        transfer_from_left(X, A)
    end
    (R_val,), (R_vec,), R_inf = eigsolve(randn(T, d, d)) do X
        transfer_from_right(X, A)
    end
    L_vec .*= sign(sum(L_vec))
    R_vec .*= sign(sum(R_vec))
    sqrtL, sqrtR = T.(sqrt(L_vec)), T.(sqrt(R_vec))
    AL = mul_matrix_from_left(mul_matrix_from_right(A, inv(sqrtL)), sqrtL)
    AR = mul_matrix_from_right(mul_matrix_from_left(A, inv(sqrtR)), sqrtR)
    C = sqrtL * sqrtR

    C ./= norm(C)
    AL ./= norm(mul_matrix_from_right(AL, C))
    AR ./= norm(mul_matrix_from_left(AR, C))
    AC = mul_matrix_from_left(AR, C)

    (AL, AR, AC, C)
end


text:
function parse_polymer(lines)
    template = popfirst!(lines)
    rules = Dict{String, Char}()
    for line in lines
        isempty(line) && continue
        pair, insertion = split(line, " -> ")
        rules[pair] = only(insertion)
    end
    polymer = Dict((pair, 0) for pair in keys(rules))
    for i in 1:length(template) - 1
        polymer[template[i:i + 1]] += 1
    end
    polymer, rules, last(template)
end

function step!(polymer, rules)
    new_polymer = Dict((pair, 0) for pair in keys(polymer))
    for (pair, n) in polymer
        insertion = rules[pair]
        new_polymer[pair[1] * insertion] += n
        new_polymer[insertion * pair[2]] += n
    end
    merge!(polymer, new_polymer)
end

function count_elements(polymer, last_element)
    element_count = Dict{Char, Int}()
    for (pair, n) in polymer
        element_count[pair[1]] = get(element_count, pair[1], 0) + n
    end
    element_count[last_element] += 1  # don't forget to count the last element in polymer
    element_count
end

polymer, rules, last_element = parse_polymer(readlines("day14/polymers.txt"))

for _ in 1:10
    step!(polymer, rules)
end

p1 = -(reverse(extrema(values(count_elements(polymer, last_element))))...)

@show p1

for _ in 1:30  # another 30 steps for a total of 40 steps
    step!(polymer, rules)
end

p2 = -(reverse(extrema(values(count_elements(polymer, last_element))))...)

@show p2


text:
# Security Considerations
# =======================
# These are tests that map directly against section 10 "Security Considerations" in the
# RFC 6455 WebSocket specification.

@testset "Security considerations" begin
    @testset "Client must choose a new masking key for each frame" begin
        # Requirement
        # @5_3-2-1 Client chooses a fresh masking key
        # @10_3-1  Client must choose a new masking key for each frame

        # Arrange
        mask = b"\x01\x02\x03\x04\x05\x06\x07\x08"
        logic, handler, writer = makeclientlogic(mask=mask)

        # Act
        text = "Hello"
        handle(logic, SendTextFrame(text, true, OPCODE_TEXT))
        handle(logic, SendTextFrame(text, true, OPCODE_TEXT))

        # Assert
        frame1 = getframe(writer, 1)
        frame2 = getframe(writer, 2)
        @test frame1.payload != frame2.payload
    end
end


text:
################################################################################
#  Copyright 2020, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################

"""
Example considering a bridge network.
"""

# load pkgs
using Unitful
using MultiStateSystems

# include the state-transition diagrams for the pipes
stdᵖ = solvedSTD(prob = [0.1, 0.9], flow = [0.0u"m^3/hr", 1.0u"m^3/hr"])

# initialize the network corresponding with the bridge network
ntwᵇⁿ = Network()

# add the sources, components and users to the network
add_source!(ntwᵇⁿ, node = 1)
add_components!(ntwᵇⁿ, edge = [(1,2), (2,4), (1,3), (3,4)],
                       std  = stdᵖ)
add_bidirectional_component!(ntwᵇⁿ, edge = (2,3), std = stdᵖ)
add_user!(ntwᵇⁿ, node = 4)

# return the network
return ntwᵇⁿ


text:
abstract type AbstractSearchResult{T} end
value(r::AbstractSearchResult) = r.value
value(r::Number) = r
location(r::AbstractSearchResult) = r.ind
payload(r::AbstractSearchResult) = r.payload
target(r::AbstractSearchResult) = r.target
targetlength(r) = lastlength(target(r))
Base.isless(a::AbstractSearchResult,b::AbstractSearchResult) = isless(value(a), value(b))
Base.isless(a::Number,b::AbstractSearchResult) = isless(value(a), value(b))
Base.isless(a::AbstractSearchResult,b::Number) = isless(value(a), value(b))
Base.promote_rule(::Type{T}, ::Type{<:AbstractSearchResult}) where T <:Number = T
Base.convert(::Type{T}, r::AbstractSearchResult) where T <:Number = T(value(r))
Base.:*(n::Number, r::AbstractSearchResult) = n*value(r)

Base.getindex(v::T, s::AbstractSearchResult{T}) where {T<:AbstractArray} =
    getwindow(v, targetlength(s), location(s))

"""
    struct BatchSearchResult{T} <: AbstractSearchResult{T}

Stores the result after applying a search over a vector of vectors. An instance `bsr` of this type acts as a number in comparisons, and vectors of BatchSearchResult can be sorted based on the `value`. `BatchSearchResult` can also act like an index into vectors. If the indexed vector `v` is of same type as `target`, then a window into `v` is returned. If the vector is of some other eltype, then `v[bsr.batch_ind]` is returned.

See accessor functions [`value`](@ref), [`location`](@ref), [`payload`](@ref), [`target`](@ref), [`targetlength`](@ref)

# Arguments:
- `target`: Contains what was searched for
- `value`: Contains the value of what's found (like the distance or cost etc.)
- `ind`: Index into the searched vector at which the `value` was found.
- `batch_ind`: Index of outer vector (or file etc.)
- `payload`: extra information associated with the search.
"""
Base.@kwdef struct BatchSearchResult{T} <: AbstractSearchResult{T}
    target::T
    value
    ind::Int
    batch_ind::Int
    payload = nothing
end

Base.:(==)(a::BatchSearchResult,b::BatchSearchResult) = a.batch_ind == b.batch_ind && location(a) == location(b)
Base.getindex(v::Vector{String}, b::BatchSearchResult) = v[b.batch_ind]


text:
# Model interconnections

"""
`series(sys1::LTISystem, sys2::LTISystem)`

Connect systems in series, equivalent to `sys2*sys1`
"""
series(sys1::LTISystem, sys2::LTISystem) = sys2*sys1

"""
`series(sys1::LTISystem, sys2::LTISystem)`

Connect systems in parallel, equivalent to `sys2+sys1`
"""
parallel(sys1::LTISystem, sys2::LTISystem) = sys1 + sys2

append() = LTISystem[]
"""
`append(systems::StateSpace...), append(systems::TransferFunction...)`

Append systems in block diagonal form
"""
function append(systems::(ST where ST<:AbstractStateSpace)...)
    ST = promote_type(typeof.(systems)...)
    Ts = systems[1].Ts
    if !all(s.Ts == Ts for s in systems)
        error("Sampling time mismatch")
    end
    A = blockdiag([s.A for s in systems]...)
    B = blockdiag([s.B for s in systems]...)
    C = blockdiag([s.C for s in systems]...)
    D = blockdiag([s.D for s in systems]...)
    return ST(A, B, C, D, Ts)
end

function append(systems::TransferFunction...)
    Ts = systems[1].Ts
    if !all(s.Ts == Ts for s in systems)
        error("Sampling time mismatch")
    end
    mat = blockdiag([s.matrix for s in systems]...)
    return TransferFunction(mat, Ts)
end

append(systems::LTISystem...) = append(promote(systems...)...)

function Base.vcat(systems::DelayLtiSystem...)
    P = vcat_1([sys.P for sys in systems]...) # See PartitionedStateSpace
    Tau = vcat([sys.Tau for sys in systems]...)
    return DelayLtiSystem(P, Tau)
end

function Base.hcat(systems::DelayLtiSystem...)
    P = hcat_1([sys.P for sys in systems]...)  # See PartitionedStateSpace
    Tau = vcat([sys.Tau for sys in systems]...)
    return DelayLtiSystem(P, Tau)
end

function Base.vcat(systems::ST...) where ST <: AbstractStateSpace
    # Perform checks
    nu = systems[1].nu
    if !all(s.nu == nu for s in systems)
        error("All systems must have same input dimension")
    end
    Ts = systems[1].Ts
    if !all(s.Ts == Ts for s in systems)
        error("Sampling time mismatch")
    end
    A = blockdiag([s.A for s in systems]...)
    B = vcat([s.B for s in systems]...)
    C = blockdiag([s.C for s in systems]...)
    D = vcat([s.D for s in systems]...)

    return ST(A, B, C, D, Ts)
end

function Base.vcat(systems::TransferFunction...)
    # Perform checks
    nu = systems[1].nu
    if !all(s.nu == nu for s in systems)
        error("All systems must have same input dimension")
    end
    Ts = systems[1].Ts
    if !all(s.Ts == Ts for s in systems)
        error("Sampling time mismatch")
    end
    mat = vcat([s.matrix for s in systems]...)
    return TransferFunction(mat, Ts)
end

Base.vcat(systems::LTISystem...) = vcat(promote(systems...)...)

function Base.hcat(systems::ST...) where ST <: AbstractStateSpace
    # Perform checks
    ny = systems[1].ny
    if !all(s.ny == ny for s in systems)
        error("All systems must have same output dimension")
    end
    Ts = systems[1].Ts
    if !all(s.Ts == Ts for s in systems)
        error("Sampling time mismatch")
    end
    A = blockdiag([s.A for s in systems]...)
    B = blockdiag([s.B for s in systems]...)
    C = hcat([s.C for s in systems]...)
    D = hcat([s.D for s in systems]...)

    return ST(A, B, C, D, Ts)
end

function Base.hcat(systems::TransferFunction...)
    # Perform checks
    ny = systems[1].ny
    if !all(s.ny == ny for s in systems)
        error("All systems must have same output dimension")
    end
    Ts = systems[1].Ts
    if !all(s.Ts == Ts for s in systems)
        error("Sampling time mismatch")
    end
    mat = hcat([s.matrix for s in systems]...)
    return TransferFunction(mat, Ts)
end

Base.hcat(systems::LTISystem...) = hcat(promote(systems...)...)

function Base._cat_t(::Val{1}, T::Type{<:LTISystem}, X...)
        vcat(convert.(T, X)...)
end

function Base._cat_t(::Val{2}, T::Type{<:LTISystem}, X...)
        hcat(convert.(T, X)...)
end

# Used in typed_hvcat
function Base.typed_hcat(::Type{T}, X...) where {T<:LTISystem}
    hcat(convert.(T, X)...)
end
# Ambiguity
Base.typed_hcat(::Type{T}, X::Number...) where {T<:LTISystem, N} = hcat(convert.(T, X)...)

# Catch special cases where inv(sys) might not be possible after promotion, like improper tf
function /(sys1::Union{StateSpace,AbstractStateSpace}, sys2::LTISystem)
    sys1new, sys2new = promote(sys1, 1/sys2)
    return sys1new*sys2new
end

# function hvcat(rows::Tuple{Vararg{Int}}, systems::Union{Number,AbstractVecOrMat{<:Number},LTISystem}...)
#     T = Base.promote_typeof(systems...)
#     nbr = length(rows)  # number of block rows
#     rs = Array{T,1}(nbr)
#     a = 1
#     for i = 1:nbr
#         rs[i] = hcat(convert.(T,systems[a:a-1+rows[i]])...)
#         a += rows[i]
#     end
#     vcat(rs...)
# end

# function _get_common_sampling_time(sys_vec::Union{AbstractVector{LTISystem},AbstractVecOrMat{<:Number},Number})
#     Ts = -1.0 # Initalize corresponding to undefined sampling time
#
#     for sys in sys_vec
#         if !all(s.Ts == Ts for s in systems])
#             error("Sampling time mismatch")
#         end
#     end
#
# end

# function Base.hcat{T<:Number}(systems::Union{T,AbstractVecOrMat{T},TransferFunction}...)
#     S = promote_type(map(e->typeof(e),systems)...) # TODO: Should be simplified
#
#     idx_first_tf = findfirst(e -> isa(e, TransferFunction), systems)
#     Ts = sys_tuple[idx_first_tf].Ts
#
#     if S <: TransferFunction
#         hcat(map(e->convert(TransferFunction,e),systems)...)
#     else
#         cat(2,systems...)
#     end
# end

blockdiag(mats::AbstractMatrix...) = blockdiag(promote(mats...)...)

function blockdiag(mats::AbstractMatrix{T}...) where T
    rows = Int[size(m, 1) for m in mats]
    cols = Int[size(m, 2) for m in mats]
    res = zeros(T, sum(rows), sum(cols))
    m = 1
    n = 1
    for ind=1:length(mats)
        mat = mats[ind]
        i = rows[ind]
        j = cols[ind]
        res[m:m + i - 1, n:n + j - 1] = mat
        m += i
        n += j
    end
    return res
end

"""
`feedback(L)` Returns L/(1+L)
`feedback(P1,P2)` Returns P1/(1+P1*P2)
"""
feedback(L::TransferFunction) = L/(1+L)
feedback(P1::TransferFunction, P2::TransferFunction) = P1/(1+P1*P2)

#Efficient implementations
function feedback(L::TransferFunction{T}) where T<:SisoRational
    if size(L) != (1,1)
        error("MIMO TransferFunction feedback isn't implemented, use L/(1+L)")
    end
    P = numpoly(L)
    Q = denpoly(L)
    tf(P, P+Q, L.Ts)
end

function feedback(L::TransferFunction{T}) where {T<:SisoZpk}
    if size(L) != (1,1)
        error("MIMO TransferFunction feedback isn't implemented, use L/(1+L)")
    end
    #Extract polynomials and create P/(P+Q)
    k = L.matrix[1].k
    denpol = numpoly(L)[1]+denpoly(L)[1]
    kden = denpol[end] # Get coeff of s^n
    # Create siso system
    sisozpk = T(L.matrix[1].z, roots(denpol), k/kden)
    return TransferFunction{T}(fill(sisozpk,1,1), L.Ts)
end

"""
`feedback(sys)`

`feedback(sys1,sys2)`

Forms the negative feedback interconnection
```julia
>-+ sys1 +-->
  |      |
 (-)sys2 +
```
If no second system is given, negative identity feedback is assumed
"""
function feedback(sys::Union{StateSpace, DelayLtiSystem})
    ninputs(sys) != noutputs(sys) && error("Use feedback(sys1, sys2) if number of inputs != outputs")
    feedback(sys,ss(Matrix{numeric_type(sys)}(I,size(sys)...)))
end

function feedback(sys1::StateSpace,sys2::StateSpace)
    if sys1.Ts != sys2.Ts # FIXME: replace with common_sample_time
        error("Sampling time mismatch")
    end

    !(iszero(sys1.D) || iszero(sys2.D)) && error("There cannot be a direct term (D) in both sys1 and sys2")
    A = [sys1.A+sys1.B*(-sys2.D)*sys1.C sys1.B*(-sys2.C);
         sys2.B*sys1.C  sys2.A+sys2.B*sys1.D*(-sys2.C)]
    B = [sys1.B; sys2.B*sys1.D]
    C = [sys1.C  sys1.D*(-sys2.C)]

    ss(A, B, C, sys1.D, sys1.Ts)
end

"""
    feedback(s1::AbstractStateSpace, s2::AbstractStateSpace;
             U1=:, Y1=:, U2=:, Y2=:, W1=:, Z1=:, W2=Int[], Z2=Int[],
             Wperm=:, Zperm=:, pos_feedback::Bool=false)

`U1`, `Y1`, `U2`, `Y2` contain the indices of the signals that should be connected.
`W1`, `Z1`, `W2`, `Z2` contain the signal indices of `s1` and `s2` that should be kept.

Specify  `Wperm` and `Zperm` to reorder [w1; w2] and [z1; z2] in the resulting statespace model.

Negative feedback is the default. Specify `pos_feedback=true` for positive feedback.

See Zhou etc. for similar (somewhat less symmetric) formulas.
"""
@views function feedback(sys1::AbstractStateSpace, sys2::AbstractStateSpace;
    U1=:, Y1=:, U2=:, Y2=:, W1=:, Z1=:, W2=Int[], Z2=Int[],
    Wperm=:, Zperm=:, pos_feedback::Bool=false)

    if sys1.Ts != sys2.Ts # FIXME: replace with common_sample_time
        error("Sampling time mismatch")
    end

    if !(isa(Y1, Colon) || allunique(Y1)); @warn "Connecting single output to multiple inputs Y1=$Y1"; end
    if !(isa(Y2, Colon) || allunique(Y2)); @warn "Connecting single output to multiple inputs Y2=$Y2"; end
    if !(isa(U1, Colon) || allunique(U1)); @warn "Connecting multiple outputs to a single input U1=$U1"; end
    if !(isa(U2, Colon) || allunique(U2)); @warn "Connecting a single output to multiple inputs U2=$U2"; end

    if (U1 isa Colon ? size(sys1, 2) : length(U1)) != (Y2 isa Colon ? size(sys2, 1) : length(Y2))
        error("Lengths of U1 ($U1) and Y2 ($Y2) must be equal")
    end
    if (U2 isa Colon ? size(sys2, 2) : length(U2)) != (Y1 isa Colon ? size(sys1, 1) : length(Y1))
        error("Lengths of U1 ($U2) and Y2 ($Y1) must be equal")
    end

    α = pos_feedback ? 1 : -1 # The sign of feedback

    s1_B1 = sys1.B[:,W1]
    s1_B2 = sys1.B[:,U1]
    s1_C1 = sys1.C[Z1,:]
    s1_C2 = sys1.C[Y1,:]
    s1_D11 = sys1.D[Z1,W1]
    s1_D12 = sys1.D[Z1,U1]
    s1_D21 = sys1.D[Y1,W1]
    s1_D22 = sys1.D[Y1,U1]

    s2_B1 = sys2.B[:,W2]
    s2_B2 = sys2.B[:,U2]
    s2_C1 = sys2.C[Z2,:]
    s2_C2 = sys2.C[Y2,:]
    s2_D11 = sys2.D[Z2,W2]
    s2_D12 = sys2.D[Z2,U2]
    s2_D21 = sys2.D[Y2,W2]
    s2_D22 = sys2.D[Y2,U2]

    if iszero(s1_D22) || iszero(s2_D22)
        A = [sys1.A + α*s1_B2*s2_D22*s1_C2        α*s1_B2*s2_C2;
                 s2_B2*s1_C2            sys2.A + α*s2_B2*s1_D22*s2_C2]

        B = [s1_B1 + α*s1_B2*s2_D22*s1_D21        α*s1_B2*s2_D21;
                      s2_B2*s1_D21            s2_B1 + α*s2_B2*s1_D22*s2_D21]
        C = [s1_C1 + α*s1_D12*s2_D22*s1_C2        α*s1_D12*s2_C2;
                      s2_D12*s1_C2           s2_C1 + α*s2_D12*s1_D22*s2_C2]
        D = [s1_D11 + α*s1_D12*s2_D22*s1_D21        α*s1_D12*s2_D21;
                      s2_D12*s1_D21           s2_D11 + α*s2_D12*s1_D22*s2_D21]
    else
        # inv seems to be better than lu
        R1 = try
            inv(α*I - s2_D22*s1_D22) # slightly faster than α*inv(I - α*s2_D22*s1_D22)
        catch
            error("Ill-posed feedback interconnection,  I - α*s2_D22*s1_D22 or I - α*s2_D22*s1_D22 not invertible")
        end

        R2 = try
            inv(I - α*s1_D22*s2_D22)
        catch
            error("Ill-posed feedback interconnection,  I - α*s2_D22*s1_D22 or I - α*s2_D22*s1_D22 not invertible")
        end

        A = [sys1.A + s1_B2*R1*s2_D22*s1_C2        s1_B2*R1*s2_C2;
                 s2_B2*R2*s1_C2            sys2.A + α*s2_B2*R2*s1_D22*s2_C2]

        B = [s1_B1 + s1_B2*R1*s2_D22*s1_D21        s1_B2*R1*s2_D21;
                     s2_B2*R2*s1_D21            s2_B1 + α*s2_B2*R2*s1_D22*s2_D21]
        C = [s1_C1 + s1_D12*R1*s2_D22*s1_C2        s1_D12*R1*s2_C2;
                     s2_D12*R2*s1_C2           s2_C1 + α*s2_D12*R2*s1_D22*s2_C2]
        D = [s1_D11 + s1_D12*R1*s2_D22*s1_D21        s1_D12*R1*s2_D21;
                     s2_D12*R2*s1_D21           s2_D11 + α*s2_D12*R2*s1_D22*s2_D21]
    end

    return StateSpace(A, B[:, Wperm], C[Zperm,:], D[Zperm, Wperm], sys1.Ts)
end

"""
`feedback2dof(P,R,S,T)` Return `BT/(AR+ST)` where B and A are the numerator and denomenator polynomials of `P` respectively
`feedback2dof(B,A,R,S,T)` Return `BT/(AR+ST)`
"""
function feedback2dof(P::TransferFunction,R,S,T)
    !issiso(P) && error("Feedback not implemented for MIMO systems")
    tf(conv(poly2vec(numpoly(P)[1]),T),zpconv(poly2vec(denpoly(P)[1]),R,poly2vec(numpoly(P)[1]),S))
end

feedback2dof(B,A,R,S,T) = tf(conv(B,T),zpconv(A,R,B,S))

"""
    lft(G, Δ, type=:l)

Lower and upper linear fractional transformation between systems `G` and `Δ`.

Specify `:l` lor lower LFT, and `:u` for upper LFT.

`G` must have more inputs and outputs than `Δ` has outputs and inputs.

For details, see Chapter 9.1 in
**Zhou, K. and JC Doyle**. Essentials of robust control, Prentice hall (NJ), 1998
"""
function lft(G, Δ, type=:l)

    if !(G.nu > Δ.ny && G.ny > Δ.nu)
        error("Must have G.nu > Δ.ny and G.ny > Δ.nu for lower/upper lft")
    end

    if type == :l
        feedback(G, Δ, U1=G.nu-Δ.ny+1:G.nu, Y1=G.ny-Δ.nu+1:G.ny, W1=1:G.ny-Δ.nu, Z1=1:G.nu-Δ.ny, pos_feedback=true)
    elseif type == :u
        feedback(G, Δ, U1=1:Δ.ny, Y1=1:Δ.nu, W1=Δ.nu+1:G.ny, Z1=Δ.nu+1:G.ny, pos_feedback=true)
    else
        error("Invalid type of lft ($type), specify type=:l (:u) for lower (upper) lft")
    end
end

"""
    starprod(sys1, sys2, dimu, dimy)

Compute the Redheffer star product.

`length(U1) = length(Y2) = dimu` and `length(Y1) = length(U2) = dimy`

For details, see Chapter 9.3 in
**Zhou, K. and JC Doyle**. Essentials of robust control, Prentice hall (NJ), 1998
"""
starprod(G1, G2, dimy::Int, dimu::Int) = feedback(G1, G2,
         U1=G1.nu-dimu+1:G1.nu, Y1=G1.ny-dimy+1:G1.ny, W1=1:G1.nu-dimu, Z1=1:G1.ny-dimy,
         U2=1:dimy, Y2=1:dimu, W2=dimy+1:G2.nu, Z2=dimu+1:G2.ny,
         pos_feedback=true)
starprod(sys1, sys2) = lft(sys1, sys2, :l)


text:
module NonstationaryPhase

using Distributions

export	PhaseParms,
		fftd,
		ifftd,
		ifftdr,
		invAℓfun,
		estϕkfun,
		unnormalized_estϕkfun,
		Cℓbiasfun,
		approxCℓbiasfun,
		Cℓvarfun,
		gen_tangentMatern,
		simNPhaseGRF,
		maternk,
		σunit_to_σpixl,
		σpixl_to_σunit,
		CZZmkfun,
		CNNkfun,
		squash,
		squash!,
		radial_power

FFTW.set_num_threads(CPU_CORES)

# This source file defines the scaled Fourier transforms used
include("gridfft.jl")

# This source file includes simple implementations, for dimension 1, of the quadratic estimate, Cvar and Cbias
# which can be used to generate estimates for non-local invariant models.
# Used in `make_figure_a.jl` and `make_figure_b.jl`
include("nonlocal_invariant_qe.jl")

##########################################################
#=
Definition of the PhaseParms Type.
Holds grid, model and planned FFT parameters for the quadratic estimate.
Allows easy argument passing.
=#
#############################################################

immutable PhaseParms{dm, T1, T2}
	# grid parameters
	period::Float64
	nside::Int64
	deltx::Float64
	deltk::Float64
	nyq::Float64
	x::Array{Array{Float64,dm},1}
	k::Array{Array{Float64,dm},1}
	r::Array{Float64,dm}
	# parameters necessary for simulation
	ξk::Array{Array{Complex{Float64},dm},1}
	ηk::Array{Array{Float64,dm},1}
	Ck:: Array{Float64,dm}
	CNNk::Array{Float64,dm}
	Cϕϕk::Array{Float64,dm}
	Eϕx::Array{Float64,dm}
	# parameters necessary for the quad est
	C1k::Array{Array{Complex{Float64},dm},1}
	C2k::Array{Array{Float64,dm},2}
	CZZmk::Array{Float64,dm}
	CZZmobsk::Array{Float64,dm} # used for weights, i.e. masked version of
	# saved plans for fast fft
	FFT::T1
	IFFT::T2            # this is the unnormalized version
	FFTconst::Float64   # these hold the normalization constants
	IFFTconst::Float64
end

"""
`PhaseParms(dm, period, nside)` constructor for PhaseParms{dm,T1,T2} type
"""
function PhaseParms(dm, period, nside)
	dm_nsides = fill(nside,dm)   # [nside,...,nside] <- dm times
	deltx     = period / nside
	deltk     = 2π / period
	nyq       = 2π / (2deltx)
	x         = [fill(NaN, dm_nsides...) for i = 1:dm]
	k         = [fill(NaN, dm_nsides...) for i = 1:dm]
	r         = fill(NaN, dm_nsides...)
	ξk        = [fill(im*NaN, dm_nsides...) for i = 1:dm]
	ηk        = [fill(NaN, dm_nsides...) for i = 1:dm]
	Ck        = fill(NaN, dm_nsides...)
	CNNk      = fill(0.0, dm_nsides...)
	Eϕx       = fill(0.0, dm_nsides...)  # default at zero
	C1k       = [fill(im*NaN, dm_nsides...) for i = 1:dm]
	C2k       = [fill(NaN, dm_nsides...) for i = 1:dm, j = 1:dm]
	CZZmk     = fill(NaN, dm_nsides...)
	CZZmobsk  = fill(NaN, dm_nsides...)
	Cϕϕk      = fill(NaN, dm_nsides...)
	FFT       = plan_fft(rand(Complex{Float64},dm_nsides...); flags = FFTW.PATIENT, timelimit = 10)
	IFFT      = plan_bfft(rand(Complex{Float64},dm_nsides...); flags = FFTW.PATIENT, timelimit = 10)
	parms     = PhaseParms{dm, typeof(FFT), typeof(IFFT)}(period, nside, deltx, deltk, nyq, x, k, r,
		ξk, ηk, Ck, CNNk, Cϕϕk, Eϕx,
		C1k, C2k, CZZmk, CZZmobsk, FFT, IFFT, (deltx / √(2π))^dm , (deltk / √(2π))^dm
		)
	parms.x[:], parms.k[:] = getgrid(parms)
	parms.r[:]  =  √(sum([abs2(kdim) for kdim in parms.k]))
	parms.ξk[:] = im.*parms.k     # default
	return parms
end
function getxkside{dm,T1,T2}(g::PhaseParms{dm,T1,T2})
	deltx    = g.period / g.nside
	deltk    = 2π / g.period
	xco_side = zeros(g.nside)
	kco_side = zeros(g.nside)
	for j in 0:(g.nside-1)
		xco_side[j+1] = (j < g.nside/2) ? (j*deltx) : (j*deltx - g.period)
		kco_side[j+1] = (j < g.nside/2) ? (j*deltk) : (j*deltk - 2*π*g.nside/g.period)
	end
	xco_side, kco_side
end
function getgrid{T1,T2}(g::PhaseParms{1,T1,T2})
	xco_side, kco_side = getxkside(g)
	xco      = Array{Float64,1}[ xco_side ]
	kco      = Array{Float64,1}[ kco_side ]
	return xco, kco
end
function meshgrid(side_x,side_y)
    	nx = length(side_x)
    	ny = length(side_y)
    	xt = repmat(vec(side_x).', ny, 1)
    	yt = repmat(vec(side_y)  , 1 , nx)
    	xt, yt
end
function getgrid{T1,T2}(g::PhaseParms{2,T1,T2})
	xco_side, kco_side = getxkside(g)
	kco1, kco2 = meshgrid(kco_side, kco_side)
	xco1, xco2 = meshgrid(xco_side, xco_side)
	kco    = Array{Float64,2}[kco1, kco2]
	xco    = Array{Float64,2}[xco1, xco2]
	return xco, kco
end

import Base.show
function Base.show{dm, T1, T2}(io::IO, parms::PhaseParms{dm, T1, T2})
	for vs in fieldnames(parms)
		(vs != :FFT) && (vs != :IFFT) && println(io, "$vs => $(getfield(parms,vs))")
		println("")
	end
end

##########################################################
#=
Closures for generating ηkfun, Ckfun, C1kfun, C2kfun
=#
#############################################################

"""
Closure which generates a nonstationay phase model in ℝ^d by specifying matern C_k and tildeC_k.
"""
function gen_tangentMatern(; ν = 2.0, ρ = 1.0, σ  = 1.0, tild_ρ = 1.0, tild_ν = 2.5, t_0 = 1.0)
	Ckfun(kco) = maternk(kco; ν=ν, ρ=ρ, σ=σ)
	function ηkfun{dm}(kco::Array{Array{Float64,dm},1})
		d0 = 4ν / ρ / ρ
		dt = 4tild_ν / tild_ρ / tild_ρ
		r2  = zero(kco[1])
		for jj = 1:dm
			r2 += abs2(kco[jj])
		end
		X0  = Distributions.Beta(dm/2, ν)
		Xt  = Distributions.Beta(dm/2, tild_ν)
		F0        = Distributions.cdf(X0, r2 ./ (r2 + d0))
		invFtF0   = (1./Distributions.quantile(Xt, F0) - 1) .^ (-1/2)
		invFtF0 .*= √(dt)
		psiprime = Array{Float64,dm}[ squash( invFtF0.*kco[jj]./√(r2) ) for jj=1:dm]
		return Array{Float64,dm}[ (psiprime[jj] .- kco[jj]) ./ t_0 for jj=1:dm]
	end
	function C1kfun{dm}(kco::Array{Array{Float64,dm},1})
		local ηkval = ηkfun(kco)
		local Ckval = Ckfun(kco)
		local  rtnkimag = Array{Complex{Float64},dm}[ im .* ηkval[j] .* Ckval ./  (2π) ^ (dm/2) for j = 1:dm]
		# Note: the  1/((2π)^(dm/2)) makes C1kfun = fft of C1x
		return rtnkimag::Array{Array{Complex{Float64},dm},1}
	end
	function C2kfun{dm}(kco::Array{Array{Float64,dm},1})
		Ck  = Ckfun(kco)
		C2k = [zero(kco[1]) for i = 1:dm, j = 1:dm]
		for ii = 1:dm, jj = 1:dm
			C2k[ii,jj]   = ηkfun(kco)[ii] .* ηkfun(kco)[jj] .* Ck
			C2k[ii,jj] ./= - 2 * (2π) ^ (dm/2)
		end
		return C2k
	end

	return ηkfun::Function, Ckfun::Function, C1kfun::Function, C2kfun::Function
end

"""
Closure which generates a nonstationay phase model in ℝ^1 by specifying a band limited matern with tangent adjustment.
"""
function gen_bandlimited_tangentMatern(nyq; ν = 2.0, ρ = 1.0, σ  = 1.0, tild_ρ = 1.0, tild_ν = 2.5, t_0 = 1.0)
	function Ckfun{dm}(kco::Array{Array{Float64,dm},1})
		r       = √( sum( Array{Float64,dm}[ abs2(kco[jj]) for jj = 1:dm ] ) )
		tanr    = ((2nyq)/π) .* tan( r .* (π/(2nyq)) )
		sec2r   = abs2( sec( r .* (π/(2nyq)) ) )
		tanargs =  Array{Float64,dm}[ squash(tanr .* kco[jj] ./ r) for jj = 1:dm]
		rtnk    = maternk(tanargs; ν=ν, ρ=ρ, σ=σ)
		rtnk  .*= sec2r
		squash!(rtnk)
		return	rtnk
	end
	function ηkfun{dm}(kco::Array{Array{Float64,dm},1})
		d0 = 4ν / ρ / ρ
		dt = 4tild_ν / tild_ρ / tild_ρ
		X0  = Distributions.Beta(dm/2, ν)
		Xt  = Distributions.Beta(dm/2, tild_ν)
		r       = √( sum( Array{Float64,dm}[ abs2(kco[jj]) for jj = 1:dm ] ) )
		tanr    = ((2nyq)/π) .* tan( r .* (π/(2nyq)) )
		sec2r   = abs2( sec( r .* (π/(2nyq)) ) )
		F0        = Distributions.cdf(X0, abs2(tanr) ./ (abs2(tanr) + d0))
		invFtF0   = (1./Distributions.quantile(Xt, F0) - 1) .^ (-1/2)
		invFtF0 .*= √(dt)
		invFtF0   = ((2nyq)/π) .* atan( invFtF0 .* (π/(2nyq)) )
		psiprime = Array{Float64,dm}[ squash(invFtF0.*kco[jj]./r) for jj=1:dm]
		return Array{Float64,dm}[ (psiprime[jj] .- kco[jj]) ./ t_0 for jj=1:dm]
	end
	function C1kfun{dm}(kco::Array{Array{Float64,dm},1})
		local ηkval = ηkfun(kco)
		local Ckval = Ckfun(kco)
		local  rtnkimag = Array{Complex{Float64},dm}[ im .* ηkval[j] .* Ckval ./  (2π) ^ (dm/2) for j = 1:dm]
		# Note: the  1/((2π)^(dm/2)) makes C1kfun = fft of C1x
		return rtnkimag::Array{Array{Complex{Float64},dm},1}
	end
	function C2kfun{dm}(kco::Array{Array{Float64,dm},1})
		Ck  = Ckfun(kco)
		C2k = [zero(kco[1]) for i = 1:dm, j = 1:dm]
		for ii = 1:dm, jj = 1:dm
			C2k[ii,jj]   = ηkfun(kco)[ii] .* ηkfun(kco)[jj] .* Ck
			C2k[ii,jj] ./= - 2 * (2π) ^ (dm/2)
		end
		return C2k
	end
	return ηkfun::Function, Ckfun::Function, C1kfun::Function, C2kfun::Function
end

##########################################################
#=
The quadratic estimate, Aℓ, Cℓvar and Cℓbias
=#
#############################################################

"""
Computes the quadratic estimate *with* normalization Aℓ
"""
function estϕkfun{dm,T1,T2}(zk, parms::PhaseParms{dm,T1,T2})
	rtnk   = unnormalized_estϕkfun(zk, parms)
    rtnk ./= invAℓfun(parms)
	squash!(rtnk)
    return rtnk
end

"""
Computes the quadratic estimate before normalization with Aℓ
"""
function unnormalized_estϕkfun{dm,T1,T2}(zk, parms::PhaseParms{dm,T1,T2})
	deltx, deltk   = parms.deltx, parms.deltk
	zkinvCZZ       = similar(zk)
	zkinvCZZ_x_2im = similar(zk)
	for inx in eachindex(zkinvCZZ)
		tmp                 = zk[inx] / parms.CZZmobsk[inx]
		zkinvCZZ[inx]       = isnan(tmp) ? Complex(0.0) : isfinite(tmp) ? tmp : Complex(0.0)
		zkinvCZZ_x_2im[inx] = isnan(tmp) ? Complex(0.0) : isfinite(tmp) ? (2*im*tmp) : Complex(0.0)
	end
	Ax   =  parms.IFFT * (zkinvCZZ)
	Bpx  = zeros(Complex{Float64}, size(zk))
	rtnk = zeros(Complex{Float64}, size(zk))
	for p = 1:dm
		Bpx[:] = parms.IFFT * (zkinvCZZ_x_2im .* imag(parms.C1k[p]))
		rtnk  += conj(parms.ξk[p]) .* (parms.FFT * (Ax .* Bpx))
	end
	scale!(rtnk, (parms.IFFTconst)^2 * parms.FFTconst)
    return rtnk
end
function unnormalized_estϕkfun{dm,T1,T2}(xk, yk, parms::PhaseParms{dm,T1,T2})
	deltx, deltk = parms.deltx, parms.deltk
	xkinvCZZ     = similar(xk)
	ykinvCZZ     = similar(yk)
	for inx in eachindex(xkinvCZZ, ykinvCZZ)
		xtmp          = xk[inx] / parms.CZZmobsk[inx]
		ytmp          = yk[inx] / parms.CZZmobsk[inx]
		xkinvCZZ[inx] = isnan(xtmp) ? Complex(0.0) : isfinite(xtmp) ? xtmp : Complex(0.0)
		ykinvCZZ[inx] = isnan(ytmp) ? Complex(0.0) : isfinite(ytmp) ? ytmp : Complex(0.0)
	end
	A1x  =  parms.IFFT * xkinvCZZ
	A2x  =  parms.IFFT * ykinvCZZ
	Dpx  = zeros(Complex{Float64}, size(xk))
	Cpx  = zeros(Complex{Float64}, size(xk))
	rtnk = zeros(Complex{Float64}, size(xk))
	for p = 1:dm
		Dpx[:] = parms.IFFT * (parms.C1k[p]  .* ykinvCZZ)
		Cpx[:] = parms.IFFT * (conj(parms.C1k[p]) .* xkinvCZZ)
		rtnk   +=  conj(parms.ξk[p]) .* (parms.FFT * (A1x .* Dpx - A2x .* Cpx))
	end
	scale!(rtnk, (parms.IFFTconst)^2 * parms.FFTconst)
    return rtnk
end

"""
Computes 1/Aℓ when CXXk == CZZmobsk
"""
function invAℓfun{dm,T1,T2}(CXXk, parms::PhaseParms{dm,T1,T2})
	deltx, deltk = parms.deltx, parms.deltk
	CXXinvCZZ = CXXk ./ parms.CZZmobsk ./ parms.CZZmobsk
	squash!(CXXinvCZZ)
	ABCDx = zeros(Complex{Float64}, size(CXXk))
	rtnk  = zeros(Complex{Float64}, size(CXXk))
	for p = 1:dm   # diag terms
		Apqx = ifftd( parms.C1k[p] .* conj(parms.C1k[p]) .* CXXinvCZZ, deltk)
		Bx   = ifftd(                                       CXXinvCZZ, deltk)
		Cpx  = ifftd(                       parms.C1k[p] .* CXXinvCZZ, deltk)
		Dpx  = ifftd(                 conj(parms.C1k[p]) .* CXXinvCZZ, deltk)
		for ix in eachindex(ABCDx)
			ABCDx[ix] = 2 * Apqx[ix] * Bx[ix] - Cpx[ix] * Cpx[ix] - Dpx[ix] * Dpx[ix]
		end
		rtnk += abs2(parms.ξk[p]) .* fftd(ABCDx, deltx)
	end
	for p = 1:dm, q = (p+1):(dm)   # off diag terms
		Apqx = ifftd( parms.C1k[p] .* conj(parms.C1k[q]) .* CXXinvCZZ, deltk)
		Bx   = ifftd(                                       CXXinvCZZ, deltk)
		Cpx  = ifftd(                       parms.C1k[p] .* CXXinvCZZ, deltk)
		Cqx  = ifftd(                       parms.C1k[q] .* CXXinvCZZ, deltk)
		Dpx  = ifftd(                 conj(parms.C1k[p]) .* CXXinvCZZ, deltk)
		Dqx  = ifftd(                 conj(parms.C1k[q]) .* CXXinvCZZ, deltk)
		for ix in eachindex(ABCDx)
			ABCDx[ix] = 2 * Apqx[ix] * Bx[ix] - Cpx[ix] * Cqx[ix] - Dpx[ix] * Dqx[ix]
		end
		rtnk += 2 * parms.ξk[p] .* conj(parms.ξk[q]) .* fftd(ABCDx, deltx)
	end
	return rtnk::Array{Complex{Float64},dm}
end
function invAℓfun{dm,T1,T2}(parms::PhaseParms{dm,T1,T2})
	CXXk = parms.CZZmk + parms.CNNk
	return invAℓfun(CXXk, parms)
end

"""
Computes Cℓvar, i.e. the variance spectral density of the quadratic estimate
"""
function Cℓvarfun{dm,T1,T2}(CXXk, parms::PhaseParms{dm,T1,T2})
	invAℓ = invAℓfun(parms)
	rtnk  = invAℓfun(CXXk, parms) ./ invAℓ ./ invAℓ
	squash!(rtnk)
	scale!(rtnk, 2 * (2π) ^ (-dm/2) )
	return real(rtnk)
end
function Cℓvarfun{dm,T1,T2}(parms::PhaseParms{dm,T1,T2})
	invAℓ = invAℓfun(parms)
	rtnk  = 1.0 ./ invAℓ
	squash!(rtnk)
	scale!(rtnk,  2 * (2π) ^ (-dm/2) )
	return real(rtnk)
end

"""
Cℓbiasfun. Function for computing the second order bias in the quadratic estimate.
"""
function Cℓbiasfun{dm,T1,T2}(parms::PhaseParms{dm,T1,T2})
	m1s   = ones(parms.Cϕϕk)
	boxbool = trues(size(parms.Cϕϕk))
	find_boxbool = find(boxbool)
    rtnkdiag  = zeros(parms.x[1])
    for p = 1:dm, q = 1:dm, p′ = 1:dm, q′ = 1:dm
    	ϕC2kℓ_pq   =  unnormalized_estϕkfun(parms.C2k[p,q],   m1s, parms)
    	ϕC2kℓ_p′q′ =  unnormalized_estϕkfun(parms.C2k[p′,q′], m1s, parms)
    	biasupdate!(p,q,p′,q′, parms,
    		ϕC2kℓ_pq,
    		ϕC2kℓ_p′q′,
    		m1s,
    		find_boxbool,
    		rtnkdiag
    		)
    end
    constant = 4 *  (parms.deltk / 2π) ^ dm
    rtnk = constant .* rtnkdiag ./ abs2(invAℓfun(parms))
    squash!(rtnk)
    return rtnk
end
function biasupdate!{dm,T1,T2}(p,q,p′,q′, parms::PhaseParms{dm,T1,T2}, ϕC2kℓ_pq, ϕC2kℓ_p′q′, m1s, find_boxbool, storage::Array{Float64,dm})
    	for ωinx in find_boxbool
    		ω = Float64[parms.k[jj][ωinx] for jj = 1:dm]
        	C2kminusω_pq      = shiftfk_by_ω(parms.C2k[p,q],   -ω, parms)
        	C2kminusω_p′q′    = shiftfk_by_ω(parms.C2k[p′,q′], -ω, parms)
    		ϕC2kℓminusω_pq    =  unnormalized_estϕkfun(C2kminusω_pq,  m1s, parms)
    		ϕC2kℓminusω_p′q′  =  unnormalized_estϕkfun(C2kminusω_p′q′, m1s, parms)
        	Cθθℓ_ωqp′   =  shiftfk_by_ω( parms.ξk[q] .* conj(parms.ξk[p′]) .* parms.Cϕϕk, -ω, parms)
        	Cθθℓ_ωqq′   =  shiftfk_by_ω( parms.ξk[q] .* conj(parms.ξk[q′]) .* parms.Cϕϕk, -ω, parms)
        	Cθθωpq′     = parms.ξk[p][ωinx] * conj(parms.ξk[q′][ωinx]) * parms.Cϕϕk[ωinx]
        	Cθθωpp′     = parms.ξk[p][ωinx] * conj(parms.ξk[p′][ωinx]) * parms.Cϕϕk[ωinx]
        	myupdateforbias!(
        		conj(ϕC2kℓ_p′q′ - ϕC2kℓminusω_p′q′),
        		ϕC2kℓ_pq - ϕC2kℓminusω_pq,
        		Cθθℓ_ωqq′,
				Cθθℓ_ωqp′,
        		Cθθωpp′,
        		Cθθωpq′,
        		storage
        		)
    	end
end
function myupdateforbias!{dm}(𝓒pq, 𝓒p′q′, Cθℓωqq′, Cθℓωqp′, Cθωpp′::Number, Cθωpq′::Number, storage::Array{Float64,dm})
	@inbounds for ind in eachindex(𝓒pq, 𝓒p′q′, Cθℓωqq′, Cθℓωqp′, storage)
		storage[ind] += real(𝓒pq[ind]  * 𝓒p′q′[ind]  * (Cθωpp′ * Cθℓωqq′[ind] +  Cθωpq′ * Cθℓωqp′[ind]))
	end
end
function shiftfk_by_ω{dm,T1,T2}(fk, ω, parms::PhaseParms{dm,T1,T2})
	fx    = parms.IFFT * fk
	ωdx  = sum([ω[jj] .* parms.x[jj] for jj = 1:dm ])
	rtnk  = parms.FFT * ( exp(-im .* ωdx) .* fx )
	scale!(rtnk, parms.IFFTconst * parms.FFTconst)
 	return rtnk
end

"""
approxCℓbiasfun
"""
function approxCℓbiasfun{dm,T1,T2}(parms::PhaseParms{dm,T1,T2})
	m1s    = ones(parms.Cϕϕk)
	∂kC2k  = make∂kC2k(parms)
	∂∂kC2k = make∂∂kC2k(parms)
	Cθθk   = makeCθθk(parms)
	rtnk   = zeros(Complex{Float64}, size(parms.k[1]))
    for p = 1:dm, q = 1:dm, p′ = 1:dm, q′ = 1:dm, ∂1 = 1:dm, ∂2 = 1:dm
    	htmp1 =      unnormalized_estϕkfun(∂kC2k[p,q,∂1]  , m1s, parms)
    	htmp2 = conj(unnormalized_estϕkfun(∂kC2k[p′,q′,∂2], m1s, parms))
    	θtmp1 = convtrans(parms.k[∂1].*parms.k[∂2].*Cθθk[p,p′], Cθθk[q,q′], parms)
    	θtmp2 = convtrans(parms.k[∂1].*parms.k[∂2].*Cθθk[p,q′], Cθθk[q,p′], parms)
    	rtnk += htmp1 .* htmp2 .* (θtmp1 + θtmp2)
    end
    for p = 1:dm, q = 1:dm, p′ = 1:dm, q′ = 1:dm, ∂1 = 1:dm, ∂2 = 1:dm, ∂3 = 1:dm
    	htmp1 =      unnormalized_estϕkfun(     ∂kC2k[p,q,∂1],       m1s, parms)
    	htmp2 = conj(unnormalized_estϕkfun(-0.5*∂∂kC2k[p′,q′,∂2,∂3], m1s, parms))
    	θtmp1 = convtrans(parms.k[∂1].*parms.k[∂2].*parms.k[∂3].*Cθθk[p,p′], Cθθk[q,q′], parms)
    	θtmp2 = convtrans(parms.k[∂1].*parms.k[∂2].*parms.k[∂3].*Cθθk[p,q′], Cθθk[q,p′], parms)
    	rtnk += htmp1 .* htmp2 .* (θtmp1 + θtmp2)
    end
    for p = 1:dm, q = 1:dm, p′ = 1:dm, q′ = 1:dm, ∂1 = 1:dm, ∂2 = 1:dm, ∂3 = 1:dm
    	htmp1 =      unnormalized_estϕkfun(-0.5*∂∂kC2k[p,q,∂1,∂2], m1s, parms)
    	htmp2 = conj(unnormalized_estϕkfun(     ∂kC2k[p′,q′,∂3]  , m1s, parms))
    	θtmp1 = convtrans(parms.k[∂1].*parms.k[∂2].*parms.k[∂3].*Cθθk[p,p′], Cθθk[q,q′], parms)
    	θtmp2 = convtrans(parms.k[∂1].*parms.k[∂2].*parms.k[∂3].*Cθθk[p,q′], Cθθk[q,p′], parms)
    	rtnk += htmp1 .* htmp2 .* (θtmp1 + θtmp2)
    end
    for p = 1:dm, q = 1:dm, p′ = 1:dm, q′ = 1:dm, ∂1 = 1:dm, ∂2 = 1:dm, ∂3 = 1:dm, ∂4 = 1:dm
    	htmp1 =      unnormalized_estϕkfun(-0.5*∂∂kC2k[p,q,∂1,∂2]   , m1s, parms)
    	htmp2 = conj(unnormalized_estϕkfun(-0.5*∂∂kC2k[p′,q′,∂3,∂4] , m1s, parms))
    	θtmp1 = convtrans(parms.k[∂1].*parms.k[∂2].*parms.k[∂3].*parms.k[∂4].*Cθθk[p,p′], Cθθk[q,q′], parms)
    	θtmp2 = convtrans(parms.k[∂1].*parms.k[∂2].*parms.k[∂3].*parms.k[∂4].*Cθθk[p,q′], Cθθk[q,p′], parms)
    	rtnk += htmp1 .* htmp2 .* (θtmp1 + θtmp2)
    end
    rtnk .*= 4 * ((2π) ^ (-dm/2))
    rtnk ./= abs2(invAℓfun(parms))
    squash!(rtnk)
    return real(rtnk)
end
# --- these are helper functions for approxCℓbias
function make∂kC2k{dm,T1,T2}(parms::PhaseParms{dm,T1,T2})
	rtnk = [zeros(Complex{Float64},size(parms.k[1])) for p = 1:dm, q = 1:dm, deriv = 1:dm]
	for p = 1:dm, q = 1:dm, deriv = 1:dm
		C2xpq  = ifftd(parms.C2k[p,q], parms.deltk)
		rtnk[p,q,deriv][:] =  fftd(-im .* parms.x[deriv] .* C2xpq, parms.deltx)
	end
	return rtnk
end
function make∂∂kC2k{dm,T1,T2}(parms::PhaseParms{dm,T1,T2})
	rtnk = [zeros(Complex{Float64},size(parms.k[1])) for p = 1:dm, q = 1:dm, ∂1 = 1:dm, ∂2 = 1:dm]
	for p = 1:dm, q = 1:dm, ∂1 = 1:dm, ∂2 = 1:dm
		C2xpq  = ifftd(parms.C2k[p,q], parms.deltk)
		rtnk[p,q,∂1,∂2][:] =  fftd( - parms.x[∂1] .* parms.x[∂2] .* C2xpq, parms.deltx)
	end
	return rtnk
end
function makeCθθk{dm,T1,T2}(parms::PhaseParms{dm,T1,T2})
	rtnk = [zeros(Complex{Float64},size(parms.k[1])) for p = 1:dm, q = 1:dm]
	for p = 1:dm, q = 1:dm
		rtnk[p,q][:] =  parms.ξk[p] .* conj(parms.ξk[q]) .* parms.Cϕϕk
	end
	return rtnk
end
function convtrans{dm,T1,T2}(Ak, Bk, parms::PhaseParms{dm,T1,T2})
	rtnk   = parms.FFT * ( (parms.IFFT * Ak) .* (parms.IFFT * Bk) )
 	rtnk .*= parms.FFTconst .* parms.IFFTconst ^ 2
 	return rtnk
end

##########################################################
#=
The marginal spectral density, CZZmkfun
=#
#############################################################

"""
Computes the marginal spectral density
"""
function CZZmkfun{dm}(Cϕϕk::Array{Float64,dm}, Ck, ηk, ξk, x, k, deltx, deltk)
	abs2k = sum([abs2(kdim) for kdim in k])
	index_xeq0 = findmin(abs2k)[2]
	Σθ = [zero(x[1]) for i = 1:dm, j = 1:dm]
	for ii = 1:dm, jj = 1:dm
		Cθix_θjx   = ifftdr( ξk[ii] .* conj(ξk[jj]) .* Cϕϕk , deltk)
		Cθix_θjx .*= (2π) ^ (-dm/2)
		Σθx_ij     = Cθix_θjx[index_xeq0] .- Cθix_θjx
		Σθx_ij   .*= 2
		Σθ[ii,jj]  = copy(Σθx_ij)
	end
	imx    = im .* x
	CZmx   = zeros(Complex{Float64}, size(Cϕϕk))
	tmp    = zeros(Complex{Float64}, size(Cϕϕk))
	tmp1   = zeros(Float64, size(Cϕϕk))
	for yi in eachindex(CZmx)
		tmp[:] = Complex(0.0)
		tmp1[:] = 0.0
		for dims1 = 1:dm
			BLAS.axpy!(imx[dims1][yi], k[dims1], tmp)
			for dims2 = 1:dm
				myscaleadd!(-0.5*Σθ[dims1,dims2][yi], ηk[dims1], ηk[dims2], tmp1)
			end
			BLAS.axpy!(1.0, tmp1, tmp)
		end
		CZmx[yi] = fastsumXexpY(Ck,tmp)
	end
	scale!((deltk ^ dm) * ((2π) ^ (-dm)), CZmx)
	CZZmk   = fftd(CZmx, deltx)
	scale!((2π) ^ (dm/2) , CZZmk)
	return abs(real(CZZmk))
end
function myscaleadd!(number, mat1, mat2, storage)
	@inbounds for ind in eachindex(mat1, mat2, storage)
		storage[ind] = storage[ind] + number * mat1[ind] * mat2[ind]
	end
end

##########################################################
#=
Simulation
=#
#############################################################
"""
# Simulate nonstationary phase model.
"""
function simNPhaseGRF{dm,T1,T2}(ϕx, parms::PhaseParms{dm,T1,T2}, parmsHR::PhaseParms{dm,T1,T2})
	ϕk     = fftd(ϕx, parms.deltx)
	boldϕx = Array{Float64,dm}[ ifftdr(parms.ξk[j] .* ϕk, parms.deltk) for j = 1:dm ]
	imboldϕx = im .* boldϕx
	imx      = im .* parms.x
	hzx_noϕ  = grfsimx(parmsHR.Ck, parmsHR.deltx, parmsHR.deltk)
	hzk_noϕ  = fftd(hzx_noϕ, parmsHR.deltx)
	zx  = zeros(Complex{Float64}, size(ϕx))
	tmp = zeros(Complex{Float64}, size(hzk_noϕ))
	for yi in eachindex(zx)
		tmp[:] = Complex(0.0)
		for dims = 1:dm   # can you loop unroll this? would it help?
			BLAS.axpy!(imboldϕx[dims][yi], parmsHR.ηk[dims], tmp)
			BLAS.axpy!(imx[dims][yi], parmsHR.k[dims], tmp)
		end
		zx[yi] = fastsumXexpY(hzk_noϕ, tmp)
	end
	zx   .*= (parmsHR.deltk ^ dm) / (2π) ^ (dm/2)
	zk     = fftd(real(zx), parms.deltx)
	zx_noϕ = downsample(hzx_noϕ, Int64(parmsHR.nside/parms.nside))
	zkobs  = zk + fftd(grfsimx(parms.CNNk, parms.deltx, parms.deltk), parms.deltx)
	return zkobs, zk, real(zx), zx_noϕ, ϕk, ϕx
end
function simNPhaseGRF{dm,T1,T2}(parms::PhaseParms{dm,T1,T2}, parmsHR::PhaseParms{dm,T1,T2})
	ϕx   = grfsimx(parms.Cϕϕk, parms.deltx, parms.deltk)
	BLAS.axpy!(1.0, parms.Eϕx, ϕx)
	zkobs, zk, zx, zx_noϕ, ϕk, ϕx = simNPhaseGRF(ϕx, parms, parmsHR)
	return zkobs, zk, zx, zx_noϕ, ϕk, ϕx
end
function fastsumXexpY(hzk, tmp) # used in simNPhaseGRF
	rtnk = Complex(0.0)
	@inbounds for ind in eachindex(hzk,tmp)
		rtnk += hzk[ind] * exp(tmp[ind])
	end
	rtnk
end
downsample{T}(mat::Array{T,1}, hfcr::Int64) =  mat[1:hfcr:end]
downsample{T}(mat::Array{T,2}, hfcr::Int64) =  mat[1:hfcr:end, 1:hfcr:end]
downsample{T}(mat::Array{T,3}, hfcr::Int64) =  mat[1:hfcr:end, 1:hfcr:end, 1:hfcr:end]
# Simulate a mean zero Gaussian random field in the pixel domain given a spectral density.
function grfsimx{T,dm}(Ckvec::Array{T,dm}, deltx, deltk)
	nsz = size(Ckvec)
	dx  = deltx ^ dm
	zzk = √(Ckvec) .* fftd(randn(nsz)./√(dx), deltx)
	return ifftdr(zzk, deltk)::Array{Float64,dm}
end

#  converting from pixel noise std to noise per-unit pixel
σunit_to_σpixl(σunit, deltx, dm) = σunit / √(deltx ^ dm)
σpixl_to_σunit(σpixl, deltx, dm) = σpixl * √(deltx ^ dm)
function CNNkfun{dm}(k::Array{Array{Float64,dm},1}, deltx; σpixl=0.0, beamFWHM=0.0)
	local absk2  = mapreduce(abs2, +, k)::Array{Float64,dm}
	local beamSQ = exp(- (beamFWHM ^ 2) * (absk2 .^ 2) ./ (8 * log(2)) )
	return ones(size(k[1])) .* σpixl_to_σunit(σpixl, deltx, dm) .^ 2 ./ beamSQ
end

function maternk{dm}(kco::Array{Array{Float64,dm},1}; ν=1.1, ρ=1.0, σ=1.0)
    d1 = 4ν / ρ / ρ
    cu = ((2π) ^ dm) * (σ ^ 2) * gamma(ν + dm/2) * ((4ν) ^ ν)
	ed = (π ^ (dm/2)) * gamma(ν) * (ρ ^ (2ν))
	# note...the extra ((2π) ^ (dm)) is so the integral equals σ^2 = ∫ C_k dk/((2π) ^ (dm))
	absk2  = mapreduce(abs2, +, kco)::Array{Float64,dm}
	rtn = (cu / ed) ./ ((d1 +  absk2) .^ (ν + dm/2))
    return rtn
end

##########################################################
#=
Cut locus
=#
#############################################################
"""
# two sided cut locus code
"""
function cutloci(parms; tmax = 20, deltatmax = 0.01)
	cutlociPos = Inf
	for t = 0:deltatmax:tmax
		if isa_cutloci(t, parms)
			cutlociPos = t
			break
		end
	end
	cutlociNeg = -Inf
	for t = 0:deltatmax:tmax
		if isa_cutloci(-t, parms)
			cutlociNeg = -t
			break
		end
	end
	return min(-cutlociNeg, cutlociPos)
end
function isa_cutloci(t, parms)
	imk   = fftshift(parms.k[1] + t*parms.ηk[1])
	return any(diff(imk) .<= 0) ? true : false
end

# this is currently un-used in but is instructive to compare with cutloci(parms)
function alternative_cutloci(parms; tmax = 20, deltatmax = 0.01)
	Ckprob   = fftshift(parms.Ck)
	Ckprob ./= sum(Ckprob)
	cutloci = Inf
	for t = 0:deltatmax:tmax
		Finv = quantile(fftshift(parms.k[1]), Ckprob)
		Fptinv = quantile(fftshift(parms.k[1] + t*parms.ηk[1]), Ckprob)
		Fntinv = quantile(fftshift(parms.k[1] - t*parms.ηk[1]), Ckprob)
		testratio = (Finv - Fptinv) ./ (Fntinv -  Fptinv)
		max_test = maximum(testratio)
		min_test = minimum(testratio)
		if (abs(max_test-min_test) > 1e-10) | (min_test < 0.) | (max_test > 1.)
			cutloci = t
			break
		end
	end
	return cutloci
end

##########################################################
#=
Miscellaneous functions
=#
#############################################################
function radial_power{dm,T1,T2}(fk, smooth::Number, parms::PhaseParms{dm,T1,T2})
	rtnk = Float64[]
	dk = parms.deltk
	kbins = collect((smooth*dk):(smooth*dk):(parms.nyq))
	for wavenumber in kbins
		indx = (wavenumber-smooth*dk) .< parms.r .<= (wavenumber+smooth*dk)
		push!(rtnk, sum(fk[indx]) / sum(indx))
	end
	return kbins, rtnk
end

squash{T<:Number}(x::T)         = isnan(x) ? zero(T) : isfinite(x) ? x : zero(T)
squash{T<:AbstractArray}(x::T)  = map(squash, x)::T
squash!{T<:AbstractArray}(x::T) = map!(squash, x)::T

function isgoodfreq{dm}(prop_of_nyq, deltx, kco::Array{Array{Float64,dm},1})
	# for generating a boolean mask
	magk = √(mapreduce(abs2, +, kco))::Array{Float64,dm}
	isgood = (magk .<= (prop_of_nyq * π / deltx)) & (magk .> 0.0)
	return isgood
end

function tent(x::Real, lend=π/2, uend=3π/2)
	# looks like this _/\_ with derivative ±1.`
	midd = lend + (uend - lend) / 2
	rtn  = 	(x ≤ lend) ? 0.0 :
			(x ≤ midd) ? ( x - lend) :
			(x ≤ uend) ? (-x + uend) : 0.0
	return rtn
end
function tent{T<:Real}(x::Array{T,1}, lend=π/2, uend=3π/2)
	return map(xr->tent(xr, lend, uend), x)
end

function Dθx_2_ϕx{dm,T1,T2}(Dθx, parms::PhaseParms{dm,T1,T2})
	Dθk = fftd(Dθx, parms.deltx)
	θk  = squash(Dθk ./ (im .* parms.k[1]))
	ϕk  = squash(θk ./ parms.ξk[1])
	ϕx = ifftdr(ϕk, parms.deltk)
	return ϕx
end
function ϕx_2_Dθx{dm,T1,T2}(ϕx, parms::PhaseParms{dm,T1,T2})
	ϕk = fftd(ϕx, parms.deltx)
	θk  = ϕk .* parms.ξk[1]
	Dθk  = θk .* (im .* parms.k[1])
	Dθx = ifftdr(Dθk, parms.deltk)
	return Dθx
end
function ϕk_2_Dθx{dm,T1,T2}(ϕk, parms::PhaseParms{dm,T1,T2})
	θk  = ϕk .* parms.ξk[1]
	Dθk  = θk .* (im .* parms.k[1])
	Dθx = ifftdr(Dθk, parms.deltk)
	return Dθx
end
function ϕk_2_θx{dm,T1,T2}(ϕk, parms::PhaseParms{dm,T1,T2})
	θk  = ϕk .* parms.ξk[1]
	θx = ifftdr(θk, parms.deltk)
	return θx
end

end # module


text:
mutable struct Pointer
    lstack::Array{UInt64,1}
    loc::Array{Int,1}
    dir::Int
    wait::Int
    printstate::Bool
    alive::Bool
    Pointer(lstack,loc,dir,wait,printstate,alive) = new(lstack,loc,dir,wait,printstate,alive)
end

mutable struct Honeycomb
    a::Array{Char,2}
    Honeycomb(a)=new(a)
end

mutable struct Time
    t::Int
    Time(t)=new(t)
end

mutable struct Debugstate
    d::Int
    Debugstate(d)=new(d)
end

"""
`beeswax(name::AbstractString)`

        beeswax("buzz.bswx") runs the program buzz.bswx up to a default limit of 1 million ticks.
        Equivalent to beeswax("buzz.bswx",0,0.0,1000000).

"""
function beeswax(name::AbstractString)
    beeswax(name,0,0.0,1000000)
end

"""
`beeswax(name::AbstractString,limit::Int)`

        beeswax("buzz.bswx",1e6) runs the program buzz.bswx up to a maximum of 1 million ticks.
        Equivalent to beeswax("buzz.bswx",0,0.0,1e6).

"""
function beeswax(name::AbstractString,limit::Int)
    beeswax(name,0,0.0,limit)
end

"""
`beeswax(name::AbstractString,debug::Debugstate,limit::Int)`

        beeswax("buzz.bswx",2,200) runs buzz.bswx in debug mode 2 up to a maximum of 200 ticks.
        Equivalent to beeswax("buzz.bswx",2,0.0,200).

    Debug modes

 *   0     no debug output.
 *   1     console output of the local stacks of all instruction pointers for every tick.
 *   2     console output of the program, instrution pointer locations and local stack contents.

"""
function beeswax(name::AbstractString,debug::Int,limit::Int)
    beeswax(name,debug,0.0,limit)
end

"""
`beeswax(name::AbstractString,debug::Debugstate,pause::Int,limit::Int)`

        beeswax("buzz.bswx",2,2,200) runs buzz.bswx in debug mode 2 up to a maximum of 200 ticks.
        Every 2 seconds the next step is displayed.
        Equivalent to beeswax("buzz.bswx",2,2.0,200).

    Debug modes

*    0     no debug output.
*    1     console output of the local stacks of all instruction pointers for every tick.
*    2     console output of the program, instrution pointer locations and local stack contents.

"""
function beeswax(name::AbstractString,debug::Int,pause::Int,limit::Int)
    beeswax(name,debug,Float64(pause),limit)
end

"""
`beeswax(name::AbstractString,debug::Debugstate,pause::Float64,limit::Int)`

        beeswax("buzz.bswx",2,0.5,200) runs buzz.bswx in debug mode 2 up to a maximum of 200 ticks.
        New steps are displayed every 0.5 seconds.

    Debug modes

 *   0     no debug output.
 *   1     console output of the local stacks of all instruction pointers for every tick.
 *   2     console output of the program, instrution pointer locations and local stack contents.

"""
function beeswax(name::AbstractString,debug::Int,pause::Float64,limit::Float64)
    beeswax(name,debug,pause,Int(limit))
end

function beeswax(name::AbstractString,debug::Int,pause::Int,limit::Float64)
    beeswax(name,debug,Float64(pause),Int(limit))
end

function beeswax(name::AbstractString,debug::Int,pause::Float64,limit::Int)

    # generate IP list
    # generate global stack
    list=Pointer[]
    sizehint!(list,10000)
    gstack=UInt64[]
    sizehint!(gstack,10000)
    ticks=Time(0)
    deb=Debugstate(debug)

    # read program file
    deb.d>2&&deb.d<9 ? println("reading program file...") :
    deb.d>=9 ? println("one-liner test mode....") : nothing
    if deb.d<9
        prog=open(name)
        code=readlines(prog)
        close(prog)
    elseif deb.d>=9
        code=split(name,"\n")
    end

    #check if it’s a valid program
    if deb.d<9
        if match(r"[\*_\\/]",readstring(name)) != nothing
            error("No starting point found. Not a valid beeswax program.")
        end
    elseif deb.d>=9
        start=0
        for i=1:length(code)
            match(r"[\*_\\/]",code[i]) != nothing ? start+=1 : nothing
        end
        start==0 ? error("No starting point found. Not a valid beeswax program.") : nothing
        deb.d==9 ? deb=Debugstate(2) : deb=Debugstate(0)
    end

    # create arena
    deb.d>2 ? println("creating honeycomb...") : nothing
    rows=length(code)
    cols=0
    for i =1:length(code)
        code[i]=chomp(code[i])
        cols=maximum([cols,textwidth(code[i])])
    end
    for i=1:length(code)
        code[i]=rpad(code[i],cols," ")
    end
    arena=Honeycomb(reshape(fill(' ',rows*cols),rows,cols))

    @inbounds for r=1:rows
        i=1;c=1
        @inbounds while i<lastindex(code[r]) || c<=length(code[r])
            try
                arena.a[r,c]=code[r][i]
            catch
            end
            i=nextind(code[r],i)
            c+=1
        end
    end

    rows=maximum(length(arena.a[:,1]))
    cols=maximum(length(arena.a[1,:]))

    # generate IPs
    deb.d>2 ? println("generating bees...") : nothing
    @inbounds for c=1:cols
        @inbounds for r=1:rows
            arena.a[r,c] == ('*')  ? pointer_all(list,r,c) :
            arena.a[r,c] == ('\\') ? pointer_diag1(list,r,c) :
            arena.a[r,c] == ('/')  ? pointer_diag2(list,r,c) :
            arena.a[r,c] == ('_')  ? pointer_horiz(list,r,c) : nothing
        end
    end

    debugger(arena,gstack,list,deb,ticks)
    ticks.t+=1

    deb.d>2 ? println("running program...") : nothing

    while length(list)>0 && ticks.t<=limit
        cleanupip(list,rows,cols)

        @inbounds for ind=length(list):-1:1
            r=list[ind].loc[1]
            c=list[ind].loc[2]
            if arena.a[r,c]==';'
                list=Pointer[];break
            elseif arena.a[r,c]==' ' && list[ind].printstate==false
                move(list,ind)
                cleanupip(list,rows,cols)
            elseif arena.a[r,c]=='X' && list[ind].printstate==false
                pointer_spread(list,ind,r,c)
                for n=length(list):-1:length(list)-3
                    move(list,n)
                end
                cleanupip(list,rows,cols)
            elseif arena.a[r,c]=='E' && list[ind].printstate==false
                pointer_spread03(list,ind,r,c)
                move(list,ind)
                move(list,length(list))
                cleanupip(list,rows,cols)
            elseif arena.a[r,c]=='H' && list[ind].printstate==false
                pointer_spread14(list,ind,r,c)
                move(list,ind)
                move(list,length(list))
                cleanupip(list,rows,cols)
            elseif arena.a[r,c]=='W' && list[ind].printstate==false
                pointer_spread25(list,ind,r,c)
                move(list,ind)
                move(list,length(list))
                cleanupip(list,rows,cols)
            else
                instruct(list,ind,gstack,arena,r,c)
                #For 'J' the jump to new location itself counts as move
                arena.a[r,c]=='J' && list[ind].printstate==false ? nothing :
                arena.a[r,c]=='J' && list[ind].printstate==true ? move(list,ind) :
                arena.a[r,c]=='v' && list[ind].printstate==true ? move(list,ind) :
                arena.a[r,c]=='v' && list[ind].printstate==false ? nothing :
                arena.a[r,c]=='^' && list[ind].printstate==true ? move(list,ind) :
                arena.a[r,c]=='^' && list[ind].printstate==false ? nothing : move(list,ind)
            end
            cleanupip(list,rows,cols)
            rows=maximum(length(arena.a[:,1]))
            cols=maximum(length(arena.a[1,:]))
        end
        pause>0.0 ? sleep(pause) : nothing
        debugger(arena,gstack,list,deb,ticks)
        ticks.t+=1
    end
    terminate()
end

function cleanupip(list::Array{Pointer,1},rows::Int,cols::Int)
    for ind=1:length(list)
        list[ind].loc[1] < 1    ? list[ind].alive=false :
        list[ind].loc[1] > rows ? list[ind].alive=false :
        list[ind].loc[2] < 1    ? list[ind].alive=false :
        list[ind].loc[2] > cols ? list[ind].alive=false : nothing
    end
    for ind=length(list):-1:1
        list[ind].alive == false ? deleteat!(list,ind) : nothing
    end
end

function move(list::Array{Pointer,1},ind::Int)
    d=list[ind].dir
    d == 0 ? list[ind].loc=list[ind].loc+[0,1] :
    d == 1 ? list[ind].loc=list[ind].loc+[-1,0] :
    d == 2 ? list[ind].loc=list[ind].loc+[-1,-1] :
    d == 3 ? list[ind].loc=list[ind].loc+[0,-1] :
    d == 4 ? list[ind].loc=list[ind].loc+[1,0] :
    d == 5 ? list[ind].loc=list[ind].loc+[1,1] : nothing
end

function slide(list::Array{Pointer,1},ind::Int,m::Array{Int,1})
    list[ind].loc=list[ind].loc+m
end

function instruct(list::Array{Pointer,1},ind::Int,gstack::Array{UInt64,1},arena::Honeycomb,row::Int,col::Int)

    c=arena.a[row,col]

    if list[ind].printstate==false
        c == '#' ? catch_all(list,ind)                :
        c == '>' ? redir_r(list,ind)                  :
        c == 'q' ? redir_dr(list,ind)                 :
        c == 'p' ? redir_dl(list,ind)                 :
        c == '<' ? redir_l(list,ind)                  :
        c == 'b' ? redir_ul(list,ind)                 :
        c == 'd' ? redir_ur(list,ind)                 :
        c == '0' ? setnum(list,ind,arena)             :
        c == '1' ? setnum(list,ind,arena)             :
        c == '2' ? setnum(list,ind,arena)             :
        c == '3' ? setnum(list,ind,arena)             :
        c == '4' ? setnum(list,ind,arena)             :
        c == '5' ? setnum(list,ind,arena)             :
        c == '6' ? setnum(list,ind,arena)             :
        c == '7' ? setnum(list,ind,arena)             :
        c == '8' ? setnum(list,ind,arena)             :
        c == '9' ? setnum(list,ind,arena)             :
        c == 'P' ? increment(list,ind)                :
        c == 'M' ? decrement(list,ind)                :
        c == '+' ? add(list,ind)                      :
        c == '-' ? sub(list,ind)                      :
        c == '.' ? mul(list,ind)                      :
        c == ':' ? intdiv(list,ind)                   :
        c == '%' ? modulus(list,ind)                  :
        c == 'B' ? pow(list,ind)                      :
        c == '&' ? bitand(list,ind)                   :
        c == '|' ? bitor(list,ind)                    :
        c == '\$' ? bitxor(list,ind)                  :
        c == '!' ? bitnot(list,ind)                   :
        c == '(' ? bitshiftleft(list,ind)             :
        c == ')' ? bitshiftright(list,ind)            :
        c == '[' ? bitrollleft(list,ind)              :
        c == ']' ? bitrollright(list,ind)             :
        c == '\'' ? iftopzero(list,ind)               :
        c == '"' ? iftopgrtzero(list,ind)             :
        c == 'K' ? iftopeqsec(list,ind)               :
        c == 'L' ? iftopgrtsec(list,ind)              :
        c == 'Q' ? skipnext(list,ind)                 :
        c == '~' ? localflip12!(list,ind)             :
        c == '@' ? localflip13!(list,ind)             :
        c == 'F' ? localallfirst!(list,ind)           :
        c == 'z' ? localallzero!(list,ind)            :
        c == 'y' ? globalrotdown!(list,ind,gstack)    :
        c == 'h' ? globalrotup!(list,ind,gstack)      :
        c == '=' ? globalduptop!(gstack)              :
        c == '?' ? globalpop!(gstack)                 :
        c == 'C' ? goutputchar(gstack)                :
        c == 'I' ? goutputint(gstack)                 :
        c == '}' ? loutputchar(list,ind)              :
        c == '{' ? loutputint(list,ind)               :
        c == '`' ? toggleoutput(list,ind)             :
        c == 'm' ? catch_diag1(list,ind)              :
        c == 'n' ? catch_diag2(list,ind)              :
        c == 'o' ? catch_diag3(list,ind)              :
        c == 'a' ? rotate_cw(list,ind)                :
        c == 'x' ? rotate_ccw(list,ind)               :
        c == 's' ? mirror_main1(list,ind)             :
        c == 't' ? mirror_main2(list,ind)             :
        c == 'u' ? mirror_main3(list,ind)             :
        c == 'j' ? mirror_half1(list,ind)             :
        c == 'k' ? mirror_half2(list,ind)             :
        c == 'l' ? mirror_half3(list,ind)             :
        c == 'O' ? mirror_all(list,ind)               :
        c == 'v' ? wait1(list,ind)                    :
        c == '^' ? wait2(list,ind)                    :
        c == 'J' ? jumpto(list,ind)                   :
        c == 'g' ? globalgetfirst(list,ind,gstack)    :
        c == 'f' ? localfirst2global(list,ind,gstack) :
        c == 'e' ? localflush2global(list,ind,gstack) :
        c == 'U' ? globalflush2local(list,ind,gstack) :
        c == 'A' ? globalstacklen!(gstack)            :
        c == 'D' ? adropto(list,ind,arena)            :
        c == 'G' ? agetfrom(list,ind,arena)           :
        c == 'Y' ? rdropto(list,ind,arena)            :
        c == 'Z' ? rgetfrom(list,ind,arena)           :
        c == 'c' ? ginputchar(gstack)                 :
        c == 'V' ? ginputstring(gstack)               :
        c == 'i' ? ginputint(gstack)                  :
        c == ',' ? linputchar(list,ind)               :
        c == 'T' ? linputint(list,ind)                :
        c == 'N' ? newline()                          :
        c == 'r' ? readfile(list,ind,gstack)          :
        c == 'w' ? writefile(list,ind,gstack)         : nothing

    else #if list[ind].printstate==true
        if arena.a[row,col] == '`'
            toggleoutput(list,ind)
        else
            print(arena.a[row,col])
        end
    end
end


text:
using Test
using CSV
using DataFrames
using StatsBase
using BenchmarkTools

using ScottishTaxBenefitModel
using ScottishTaxBenefitModel.GeneralTaxComponents
using ScottishTaxBenefitModel.STBParameters
using ScottishTaxBenefitModel.Runner: do_one_run!
using ScottishTaxBenefitModel.RunSettings: Settings, MT_Routing
using .Utils
using .ExampleHelpers

settings = Settings()

BenchmarkTools.DEFAULT_PARAMETERS.seconds = 120
BenchmarkTools.DEFAULT_PARAMETERS.samples = 2

function basic_run( ; print_test :: Bool, mtrouting :: MT_Routing  )
    settings.means_tested_routing = mtrouting
    settings.run_name="run-$(mtrouting)-$(date_string())"
    sys = [get_system(scotland=false), get_system( scotland=true )]
    results = do_one_run!( settings, sys )
end

@testset "basic run timing" begin
    for mt in instances( MT_Routing )
        println( "starting run using $mt routing")
        @time basic_run( print_test=true, mtrouting = mt )
    end
    # @benchmark frames =
    # print(t)
end

#=
if print_test
    summary_output = summarise_results!( results=results, base_results=base_results )
    print( "   deciles = $( summary_output.deciles)\n\n" )
    print( "   poverty_line = $(summary_output.poverty_line)\n\n" )
    print( "   inequality = $(summary_output.inequality)\n\n" )
    print( "   poverty = $(summary_output.poverty)\n\n" )
    print( "   gainlose_by_sex = $(summary_output.gainlose_by_sex)\n\n" )
    print( "   gainlose_by_thing = $(summary_output.gainlose_by_thing)\n\n" )
    print( "   metr_histogram= $(summary_output.metr_histogram)\n\n")
    println( "SUMMARY OUTPUT")
    println( summary_output )
    println( "as JSON")
    println( JSON.json( summary_output ))
end

=#