text:
# returns the list of subexpressions which a given tape depends on directly
function list_subexpressions(nd::Vector{NodeData})

    # TODO: code overlap with compute_gradient_sparsity
    indices = Set{Int}()
    for k in 1:length(nd)
        nod = nd[k]
        if nod.nodetype == SUBEXPRESSION
            push!(indices, nod.index)
        end
    end

    return sort(collect(indices))
end
# Order the subexpressions which main_expressions depend on such that we can
# run forward mode in this order.
function order_subexpressions(main_expressions::Vector{Vector{NodeData}},
                              subexpressions::Vector{Vector{NodeData}})
    num_sub = length(subexpressions)
    computed = falses(num_sub)
    dependencies = Dict{Int,Vector{Int}}()
    to_visit = collect(num_sub + 1 : num_sub + length(main_expressions))
    # For each subexpression k, the indices of the main expressions that depend
    # on k, possibly transitively.
    depended_on_by = [Set{Int}() for i in 1:num_sub]

    while !isempty(to_visit)
        idx = pop!(to_visit)
        if idx > num_sub
            subexpr = list_subexpressions(main_expressions[idx - num_sub])
        else
            computed[idx] && continue
            subexpr = list_subexpressions(subexpressions[idx])
            computed[idx] = true
        end
        dependencies[idx] = subexpr
        for k in subexpr
            if idx > num_sub
                push!(depended_on_by[k], idx - num_sub)
            end
            push!(to_visit, k)
        end
    end

    # Now order dependencies.
    I = Int[]
    J = Int[]
    for (idx, subexpr) in dependencies
        for k in subexpr
            push!(I, idx)
            push!(J, k)
        end
    end
    N = num_sub + length(main_expressions)
    sp = sparse(I, J, ones(length(I)), N, N)
    cmap = Vector{Int}(undef, N)
    order = reverse(Coloring.reverse_topological_sort_by_dfs(sp.rowval,
                                                             sp.colptr, N,
                                                             cmap)[1])
    # Remove the subexpressions which never appear anywhere and the indices of
    # the main expressions.
    condition(idx) = idx <= num_sub && computed[idx]
    order_filtered = collect(filter(condition, order))

    # Propagate transitive dependencies.
    for o in Iterators.reverse(order_filtered)
        @assert !isempty(depended_on_by[o])
        for k in list_subexpressions(subexpressions[o])
            union!(depended_on_by[k], depended_on_by[o])
        end
    end

    # Generate an individual order for each main expression.
    individual_order = [Int[] for i in 1:length(main_expressions)]
    for o in order_filtered
        for i in depended_on_by[o]
            push!(individual_order[i], o)
        end
    end

    return order_filtered, individual_order
end

export list_subexpressions, order_subexpressions


text:
#####
##### MixtapeInterpreter
#####

# User-extended context allows parametrization of the pipeline through
# our subtype of AbstractInterpreter
abstract type CompilationContext end

@doc(
"""
    abstract type CompilationContext end

Parametrize the Mixtape pipeline by inheriting from `CompilationContext`. Similar to the context objects in [Cassette.jl](https://julia.mit.edu/Cassette.jl/stable/contextualdispatch.html). By using the interface methods `show_after_inference`, `show_after_optimization`, `debug`, `transform`, and `optimize!` -- the user can control different parts of the compilation pipeline.
""", CompilationContext)

show_after_inference(ctx::CompilationContext) = false

@doc(
"""
    show_after_inference(ctx::CompilationContext)::Bool

Turns on a pipeline feature which will dump out `CodeInfo` after inference (including the user-defined `transform` transformation, if applied).
""", show_after_inference)

show_after_optimization(ctx::CompilationContext) = false

@doc(
"""
    show_after_optimization(ctx::CompilationContext)::Bool

Turns on a pipeline feature which will dump out `CodeInfo` after optimization (including the user-defined `optimize!` transformation, if applied).
""", show_after_optimization)

debug(ctx::CompilationContext) = false

@doc(
"""
    debug(ctx::CompilationContext)::Bool

Turn on debug tracing for inference and optimization. Displays an instrumentation trace as inference and optimization proceeds.
""", debug)

transform(ctx::CompilationContext, b) = b

@doc(
"""
    transform(ctx::CompilationContext, b::CodeInfoTools.Builder)::CodeInfoTools.Builder

User-defined transform which can operate on lowered `CodeInfo` in the form of a `CodeInfoTools.Builder` object.

Transforms might typically follow a simple "replace" format:

```julia
function transform(::MyCtx, b)
    for (k, st) in b
        replace!(b, k, swap(st))
    end
    return b
end
```

but more advanced formats are possible. For further utilities, please see [CodeInfoTools.jl](https://github.com/femtomc/CodeInfoTools.jl).
""", transform)

optimize!(ctx::CompilationContext, ir) = ir

@doc(
"""
    optimize!(ctx::CompilationContext, ir::Core.Compiler.IRCode)::Core.Compiler.IRCode

User-defined transform which can operate on inferred `Core.Compiler.IRCode`. This transform operates after a set of optimizations which mimic Julia's pipeline.
""", optimize!)

allow(f::C, args...) where {C <: CompilationContext} = false
function allow(ctx::CompilationContext, mod::Module, fn, args...)
    return allow(ctx, mod) || allow(ctx, fn, args...)
end

@doc(
"""
    allow(f::CompilationContext, args...)::Bool

Determines whether a user-defined `transform` or `optimize!` is allowed to look at a lowered `CodeInfoTools.Builder` object or `Core.Compiler.IRCode`.

The user is allowed to greenlight modules:

```julia
allow(::MyCtx, m::Module) == m == SomeModule
```

or even specific signatures

```julia
allow(::MyCtx, fn::typeof(rand), args...) = true
```
""", allow)

macro ctx(properties, expr)
    @assert(@capture(expr, struct Name_
                         body__
                     end))
    @assert(properties.head == :tuple)
    properties = properties.args
    ex = Expr(:block,
              quote
                  import Mixtape: allow, show_after_inference, show_after_optimization,
                                  debug, transform, optimize!
              end, quote
                  struct $Name <: Mixtape.CompilationContext
                      $(body...)
                  end
              end, quote
                  show_after_inference(::$Name) = $(properties[1])
                  show_after_optimization(::$Name) = $(properties[2])
                  debug(::$Name) = $(properties[3])
              end)
    ex = postwalk(rmlines ∘ unblock, ex)
    return esc(ex)
end

@doc(
"""
    @ctx(properties, expr)

Utility macro which expands to implement a subtype of `CompilationContext`. Also allow the user to easily configure the `debug` static tracing features and `show_after_inference`/`show_after_optimization` Boolean function flags.

Usage:
```julia
@ctx (b1::Bool, b2::Bool, b3::Bool) struct MyCtx
    fields...
end
```

Expands to:
```julia
struct MyCtx <: CompilationContext
    fields...
end
show_after_inference(::MyCtx) = b1
show_after_optimization(::MyCtx) = b2
debug(::MyCtx) = b3
```
""", :(@Mixtape.ctx))

struct MixtapeInterpreter{Ctx<:CompilationContext,Inner<:AbstractInterpreter} <:
       AbstractInterpreter
    ctx::Ctx
    inner::Inner
    errors::Any
    function MixtapeInterpreter(ctx::Ctx,
                                interp::Inner) where {Ctx<:CompilationContext,Inner}
        return new{Ctx,Inner}(ctx, interp, Any[])
    end
end
Base.push!(mxi::MixtapeInterpreter, e) = push!(mxi.errors, e)

get_world_counter(mxi::MixtapeInterpreter) = get_world_counter(mxi.inner)
get_inference_cache(mxi::MixtapeInterpreter) = get_inference_cache(mxi.inner)
InferenceParams(mxi::MixtapeInterpreter) = InferenceParams(mxi.inner)
OptimizationParams(mxi::MixtapeInterpreter) = OptimizationParams(mxi.inner)
Core.Compiler.may_optimize(mxi::MixtapeInterpreter) = true
Core.Compiler.may_compress(mxi::MixtapeInterpreter) = true
Core.Compiler.may_discard_trees(mxi::MixtapeInterpreter) = true
function Core.Compiler.add_remark!(mxi::MixtapeInterpreter, sv::InferenceState, msg)
    return Core.Compiler.add_remark!(mxi.inner, sv, msg)
end
lock_mi_inference(mxi::MixtapeInterpreter, mi::MethodInstance) = nothing
unlock_mi_inference(mxi::MixtapeInterpreter, mi::MethodInstance) = nothing

ext: jl


text:
export structural, local_level, linear_trend, regression, tvm

"""
    structural(y::VecOrMat{Typ}, s::Int; X::VecOrMat{Typ} = Matrix{Typ}(undef, 0, 0)) where Typ <: Real

Build state-space system for a given structural model with observations `y`, seasonality `s`, and, optionally, exogenous variables `X`.

If `y` is provided as an `Array{Typ, 1}` it will be converted to an `Array{Typ, 2}` inside the `StateSpaceModel`. The same will happen to X,
if an `Array{Typ, 1}` it will be converted to an `Array{Typ, 2}` inside the `StateSpaceModel`.
"""
function structural(y::VecOrMat{Typ}, s::Int; X::VecOrMat{Typ} = Matrix{Typ}(undef, 0, 0)) where Typ <: Real

    # Number of observations and endogenous variables
    y = y[:, :]
    n, p = size(y)

    # Number of observations and exogenous variables
    X = X[:, :]
    n_exp, p_exp = size(X)

    if p_exp > 0 && n_exp < n
        error("Number of observations in X must be greater than or equal to the number of observations in y.")
    end

    m = (1 + s + p_exp)*p
    r = 3*p

    # Observation equation
    N = max(n, n_exp)
    if p_exp > 0 # exogenous variables: Z is time-variant
        Z = Array{Typ, 3}(undef, p, m, N)
        for t = 1:N
            Z[:, :, t] = kron(
                Matrix{Typ}(I, p, p),
                [
                    X[t, :]' 1 0 1 zeros(Typ, 1, s - 2)
                ]
                )
        end
    else # no exogenous variables: Z is time-invariant
        Z = kron(
            Matrix{Typ}(I, p, p),
            [
                1 0 1 zeros(Typ, 1, s - 2)
            ]
            )
    end

    # State equation
    if p_exp > 0
        T0 = [Matrix{Typ}(I, p_exp, p_exp) zeros(Typ, p_exp, 1 + s)]
        T = kron(
            Matrix{Typ}(I, p, p),
            [
                T0;
                zeros(Typ, 1, p_exp) 1 1 zeros(Typ, 1, s - 1);
                zeros(Typ, 1, p_exp) 0 1 zeros(Typ, 1, s - 1);
                zeros(Typ, 1, p_exp) 0 0 -ones(Typ, 1, s - 1);
                zeros(Typ, s - 2, p_exp) zeros(Typ, s - 2, 2) Matrix{Typ}(I, s - 2, s - 2) zeros(Typ, s - 2)
            ]
            )
    else
        T = kron(
            Matrix{Typ}(I, p, p),
            [
                1 1 zeros(Typ, 1, s - 1);
                0 1 zeros(Typ, 1, s - 1);
                0 0 -ones(Typ, 1, s - 1);
                zeros(Typ, s - 2, 2) Matrix{Typ}(I, s - 2, s - 2) zeros(Typ, s - 2)
            ]
            )
        end

    R = kron(
        Matrix{Typ}(I, p, p),
        [
            zeros(Typ, p_exp, 3);
            Matrix{Typ}(I, 3, 3);
            zeros(Typ, s - 2, 3)
        ]
        )

    return StateSpaceModel(y, Z, T, R)

end

"""
    local_level(y::VecOrMat{Typ}) where Typ <: Real

Build state-space system for a local level model with observations y.

If `y` is provided as an `Array{Typ, 1}` it will be converted to an `Array{Typ, 2}` inside the `StateSpaceModel`.
"""
function local_level(y::VecOrMat{Typ}) where Typ <: Real

    # Number of observations and endogenous variables
    y = y[:, :]
    n, p = size(y)

    m = p
    r = p

    # Observation equation
    Z = Matrix{Typ}(I, p, p)

    # State equation
    T = Matrix{Typ}(I, p, p)
    R = Matrix{Typ}(I, p, p)

    return StateSpaceModel(y, Z, T, R)
end

"""
    linear_trend(y::VecOrMat{Typ}) where Typ <: Real

Build state-space system for a linear trend model with observations y.

If `y` is provided as an `Array{Typ, 1}` it will be converted to an `Array{Typ, 2}` inside the `StateSpaceModel`.
"""
function linear_trend(y::VecOrMat{Typ}) where Typ <: Real

    # Number of observations and endogenous variables
    y = y[:, :]
    n, p = size(y)

    m = 2*p
    r = 2*p

    # Observation equation
    Z = kron(Matrix{Typ}(I, p, p),[1 0])

    # State equation
    T = kron(Matrix{Typ}(I, p, p),[1 1; 0 1])
    R = kron(Matrix{Typ}(I, p, p),[1 0; 0 1])

    return StateSpaceModel(y, Z, T, R)
end

"""
    regression(y::VecOrMat{Typ}, X::VecOrMat{Typ}) where Typ <: Real

Build state-space system for estimating a regression model ``y_t = X_t\\beta_t + \\varepsilon_t``. Once the model is estimated
the user can recover the parameter ``\\hat{\\beta}`` by the querying the smoothed states of the model.
"""
function regression(y::VecOrMat{Typ}, X::VecOrMat{Typ}) where Typ <: Real
    # Certify that they are matrices
    y = ensure_is_matrix(y)
    X = ensure_is_matrix(X)
    # Assert sizes are adequate
    n_y, p_y = size(y)
    n_X, p_X = size(X)
    @assert n_X == n_y

    if p_y > 1
        error("StateSpaceModels currently supports regression only for univariate cases.")
    end

    # Fill Z
    Z = Array{Typ, 3}(undef, p_y, p_X, n_X)
    for t in 1:n_X
        Z[:, :, t] = X[t, :]
    end
    T = Matrix{Typ}(I, p_X, p_X)
    R = zeros(Typ, p_X, 1)
    Q = zeros(Typ, 1, 1)
    d = zeros(Typ, n_y, p_y)
    c = zeros(Typ, n_X, p_X)
    # H is the only variable to estimate
    H = build_H(p_y, Typ)

    return StateSpaceModel(y, Z, T, R, d, c, H, Q)
end

"""
    tvm(y::VecOrMat{Typ}, X::VecOrMat{Typ}) where Typ <: Real

Build state-space system for estimating a time varying regression model ``y_t = X_t\\beta_t + \\varepsilon_t``. Once the model is estimated
the user can recover the parameter ``\\hat{\\beta}`` by the querying the smoothed states of the model.
"""

function tvm(y::VecOrMat{Typ}, X::VecOrMat{Typ}; X_var = zeros(Typ, size(X, 2))) where Typ <: Real

    #Time varying
    # Certify that they are matrices
    y = ensure_is_matrix(y)
    X = ensure_is_matrix(X)
    # Assert sizes are adequate
    n_y, p_y = size(y)
    n_X, p_X = size(X)

    #We are creating a m-dimensional state-space model
    m = p_X + p_y
    @assert n_X == n_y

    if p_y > 1
        error("StateSpaceModels currently supports regression only for univariate cases.")
    end

    # Fill Z
    Z = Array{Typ, 3}(undef, p_y, m, n_X)
    for t in 1:n_X
        Z[:, :, t] = vcat(Matrix{Typ}(I, p_y, p_y), X[t, :])
    end
    T = Matrix{Typ}(I, m, m)
    R = diagm(vcat([1.0], X_var))

    Q = diagm(vcat([1.0], X_var))
    Q[findall(isequal(1.0), Q)] .= NaN

    d = zeros(Typ, n_y, p_y)
    c = zeros(Typ, n_X, m)

    H = build_H(p_y, Typ)

    return StateSpaceModel(y, Z, T, R, d, c, H, Q)
end

max_stars_repo_name: PaulMainwood/StateSpaceCPG

max_issues_repo_name: PaulMainwood/StateSpaceCPG

max_forks_repo_name: PaulMainwood/StateSpaceCPG


text:
@testset "EGPpower" begin
    @test_throws ArgumentError EGPpower(-1, .2, 2)
    @test_throws ArgumentError EGPpower(-1, .2, -2)

    σ = 1
    ξ = .3
    κ = 2
    model = EGPpower(σ, ξ, κ)

    @test minimum(model) == 0
    @test maximum(model) == Inf

    @test insupport(model, -1) == false
    @test insupport(model, 0) == true
    @test insupport(model, 1000) == true

    @test params(model) == (σ, ξ, κ)

    @test partype(model) == Float64

    @test_throws AssertionError quantile(model, 0)
    @test_throws AssertionError quantile(model, 1)
end


text:
using LineNumbers
using Base.Test

code = "f() = 1 + 2 - 3\n"
@test SourceFile(code)[1:2] == Vector{UInt8}["f() = 1 + 2 - 3".data, "".data]
data2 = Vector{UInt8}[copy(code.data),[]]
LineBreaking(UInt64(1),SourceFile(code),data2)[7] = 'x'
@test bytestring(data2[1]) == "f() = x + 2 - 3\n"
@test length(SourceFile("f() = 1 + 2 - 3")) == 1


text:
mutable struct HSAExecutable{B} # FIXME: Don't parameterize on Mem.Buffer
    agent::HSAAgent
    executable::Ref{HSA.Executable}
    data::Vector{UInt8}
    globals::Dict{Symbol,B}
end

### @cfunction Callbacks ###

function iterate_exec_agent_syms_cb(exe::HSA.Executable, agent::HSA.Agent,
                                    sym::HSA.ExecutableSymbol,
                                    sym_ref::Ptr{Cvoid})

    sym_ref = Base.unsafe_convert(Ptr{HSA.ExecutableSymbol}, sym_ref)
    sym_type = Ref{HSA.SymbolKind}()
    getinfo(sym, HSA.EXECUTABLE_SYMBOL_INFO_TYPE,
                                   sym_type) |> check

    if sym_type[] == HSA.SYMBOL_KIND_KERNEL
        len = Ref(0)
        getinfo(sym, HSA.EXECUTABLE_SYMBOL_INFO_NAME_LENGTH, len) |> check
        name = Vector{UInt8}(undef, len[])
        getinfo(sym, HSA.EXECUTABLE_SYMBOL_INFO_NAME, name) |> check
        @debug "    Symbol Name: $(String(name))"
        Base.unsafe_store!(sym_ref, sym)
    end

    return HSA.STATUS_SUCCESS
end

# TODO Docstring
function HSAExecutable(agent::HSAAgent, data::Vector{UInt8}, symbol::String; globals=())
    #= NOTE: Everything I can see indicates that profile is always FULL
    profile = Ref{HSA.Profile}()
    HSA.agent_get_info(agent.agent, HSA.AGENT_INFO_PROFILE, profile) |> check
    =#

    code_object_reader = Ref{HSA.CodeObjectReader}(HSA.CodeObjectReader(0))
    HSA.code_object_reader_create_from_memory(data, sizeof(data),
                                              code_object_reader) |> check

    executable = Ref{HSA.Executable}()
    HSA.executable_create_alt(HSA.PROFILE_BASE,
                              HSA.DEFAULT_FLOAT_ROUNDING_MODE_NEAR,
                              C_NULL, executable) |> check

    _globals = Dict{Symbol,Mem.Buffer}()
    for (gbl,sz) in globals
        gbl_buf = Mem.alloc(agent, sz; coherent=true)
        HSA.executable_agent_global_variable_define(executable[], agent.agent,
                                                    string(gbl), gbl_buf.ptr) |> check
        _globals[gbl] = gbl_buf
    end

    HSA.executable_load_agent_code_object(executable[], agent.agent,
                                          code_object_reader[],
                                          C_NULL, C_NULL) |> check

    HSA.executable_freeze(executable[], "") |> check

    exe = HSAExecutable(agent, executable, data, _globals)

    # TODO: Ensure no derived kernels are in flight during finalization
    hsaref!()
    finalizer(exe) do exe
        HSA.executable_destroy(exe.executable[]) |> check
        for buf in values(exe.globals)
            Mem.free(buf)
        end
        HSA.code_object_reader_destroy(code_object_reader[]) |> check
        hsaunref!()
    end

    return exe
end

function get_global(exe::HSAExecutable, symbol::Symbol)
    @assert symbol in keys(exe.globals) "No such global in executable: $symbol"
    return exe.globals[symbol]::Mem.Buffer
end


text:
@testset "Dense-Sparse in-place speed bug" begin
    # https://discourse.julialang.org/t/asymmetric-speed-of-in-place-sparse-dense-matrix-product/10256/5
    A = sprand(1000,1000,0.01);
    B = rand(1000,1000);
    C = similar(B);

    # check that the order doesn't matter much
    sparse_dense = @belapsed mul!($C,$A,$B);
    dense_sparse = @belapsed mul!($C,$A,$B);
    @test abs(sparse_dense - dense_sparse)/(sparse_dense+dense_sparse) < 0.7
    # let's hope that this thresh works on other machines
end

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
using DistributionsAD

using Combinatorics
using Distributions
using PDMats

using Random, LinearAlgebra, Test

using Distributions: meanlogdet
using DistributionsAD: TuringUniform, TuringMvNormal, TuringMvLogNormal,
                       TuringPoissonBinomial, TuringDirichlet
using StatsBase: entropy
using StatsFuns: StatsFuns, logsumexp, logistic

Random.seed!(1) # Set seed that all testsets should reset to.

const GROUP = get(ENV, "GROUP", "All")

include("test_utils.jl")

if GROUP == "All" || GROUP == "Others"
    include("others.jl")
end

if GROUP == "All" || GROUP in ("ForwardDiff", "Zygote", "ReverseDiff", "Tracker")
    include("ad/utils.jl")
    include("ad/others.jl")
    include("ad/chainrules.jl")
    include("ad/distributions.jl")
end


text:
# Initializing Python (surprisingly complicated; see also deps/build.jl)

#########################################################################

# global PyObject constants, initialized to NULL and then overwritten in __init__
# (eventually, the ability to define global const in __init__ may go away,
#  and in any case this is better for type inference during precompilation)
const inspect = PyNULL()
const builtin = PyNULL()
const BuiltinFunctionType = PyNULL()
const MethodType = PyNULL()
const MethodWrapperType = PyNULL()
const ufuncType = PyNULL()
const format_traceback = PyNULL()
const pyproperty = PyNULL()
const jlfun2pyfun = PyNULL()
const c_void_p_Type = PyNULL()

# other global constants initialized at runtime are defined via Ref
# or are simply left as non-const values
const pynothing = Ref{PyPtr}(0)
const pyxrange = Ref{PyPtr}(0)

#########################################################################
# initialize jlWrapType for pytype.jl

function pyjlwrap_init()
    # PyMemberDef stores explicit pointers, hence must be initialized at runtime
    push!(pyjlwrap_members, PyMemberDef(pyjlwrap_membername,
                                        T_PYSSIZET, sizeof_PyObject_HEAD, READONLY,
                                        pyjlwrap_doc),
                            PyMemberDef(C_NULL,0,0,0,C_NULL))

    # all cfunctions must be compiled at runtime
    pyjlwrap_dealloc_ptr = @cfunction(pyjlwrap_dealloc, Cvoid, (PyPtr,))
    pyjlwrap_repr_ptr = @cfunction(pyjlwrap_repr, PyPtr, (PyPtr,))
    pyjlwrap_hash_ptr = @cfunction(pyjlwrap_hash, UInt, (PyPtr,))
    pyjlwrap_hash32_ptr = @cfunction(pyjlwrap_hash32, UInt32, (PyPtr,))
    pyjlwrap_call_ptr = @cfunction(pyjlwrap_call, PyPtr, (PyPtr,PyPtr,PyPtr))
    pyjlwrap_getattr_ptr = @cfunction(pyjlwrap_getattr, PyPtr, (PyPtr,PyPtr))
    pyjlwrap_getiter_ptr = @cfunction(pyjlwrap_getiter, PyPtr, (PyPtr,))

    # detect at runtime whether we are using Stackless Python
    try
        pyimport("stackless")
        Py_TPFLAGS_HAVE_STACKLESS_EXTENSION[] = Py_TPFLAGS_HAVE_STACKLESS_EXTENSION_
    catch
    end

    PyTypeObject!(jlWrapType, "PyCall.jlwrap", sizeof(Py_jlWrap)) do t::PyTypeObject
        t.tp_flags |= Py_TPFLAGS_BASETYPE
        t.tp_members = pointer(pyjlwrap_members);
        t.tp_dealloc = pyjlwrap_dealloc_ptr
        t.tp_repr = pyjlwrap_repr_ptr
        t.tp_call = pyjlwrap_call_ptr
        t.tp_getattro = pyjlwrap_getattr_ptr
        t.tp_iter = pyjlwrap_getiter_ptr
        t.tp_hash = sizeof(Py_hash_t) < sizeof(Int) ?
                    pyjlwrap_hash32_ptr : pyjlwrap_hash_ptr
    end
end

#########################################################################

function __init__()
    # sanity check: in Pkg for Julia 0.7+, the location of Conda can change
    # if e.g. you checkout Conda master, and we'll need to re-build PyCall
    # for something like pyimport_conda to continue working.
    if conda && dirname(python) != abspath(Conda.PYTHONDIR)
        error("Using Conda.jl python, but location of $python seems to have moved to $(Conda.PYTHONDIR).  Re-run Pkg.build(\"PyCall\") and restart Julia.")
    end

    # issue #189
    libpy_handle = libpython === nothing ? C_NULL :
        Libdl.dlopen(libpython, Libdl.RTLD_LAZY|Libdl.RTLD_DEEPBIND|Libdl.RTLD_GLOBAL)

    already_inited = 0 != ccall((@pysym :Py_IsInitialized), Cint, ())

    if !already_inited
        Py_SetPythonHome(libpy_handle, PYTHONHOME, wPYTHONHOME, pyversion)
        if !isempty(pyprogramname)
            if pyversion.major < 3
                ccall((@pysym :Py_SetProgramName), Cvoid, (Cstring,), pyprogramname)
            else
                ccall((@pysym :Py_SetProgramName), Cvoid, (Ptr{Cwchar_t},), wpyprogramname)
            end
        end
        ccall((@pysym :Py_InitializeEx), Cvoid, (Cint,), 0)
    end

    # Will get reinitialized properly on first use
    Compat.Sys.iswindows() && (PyActCtx[] = C_NULL)

    # Make sure python wasn't upgraded underneath us
    new_pyversion = vparse(split(unsafe_string(ccall(@pysym(:Py_GetVersion),
                               Ptr{UInt8}, ())))[1])

    if new_pyversion.major != pyversion.major
        error("PyCall precompiled with Python $pyversion, but now using Python $new_pyversion; ",
              "you need to relaunch Julia and re-run Pkg.build(\"PyCall\")")
    end

    copy!(inspect, pyimport("inspect"))
    copy!(builtin, pyimport(pyversion.major < 3 ? "__builtin__" : "builtins"))
    copy!(pyproperty, pybuiltin(:property))

    pyexc_initialize() # mappings from Julia Exception types to Python exceptions

    # cache Python None -- PyPtr, not PyObject, to prevent it from
    # being finalized prematurely on exit
    pynothing[] = @pyglobalobj(:_Py_NoneStruct)

    # xrange type (or range in Python 3)
    pyxrange[] = @pyglobalobj(:PyRange_Type)

    # ctypes.c_void_p for Ptr types
    copy!(c_void_p_Type, pyimport("ctypes")["c_void_p"])

    # traceback.format_tb function, for show(PyError)
    copy!(format_traceback, pyimport("traceback")["format_tb"])

    init_datetime()
    pyjlwrap_init()

    # jlwrap is a class, and when assigning it to an object
    #    obj[:foo] = some_julia_function
    # it won't behave like a regular Python method because it's not a Python
    # function (in particular, `self` won't be passed to it). The solution is:
    #    obj[:foo] = jlfun2pyfun(some_julia_function)
    # This is a bit of a kludge, obviously.
    copy!(jlfun2pyfun,
          pyeval_("""lambda f: lambda *args, **kwargs: f(*args, **kwargs)"""))

    if !already_inited
        # some modules (e.g. IPython) expect sys.argv to be set
        @static if VERSION >= v"0.7.0-DEV.1963"
            ref0 = Ref{UInt32}(0)
            GC.@preserve ref0 ccall(@pysym(:PySys_SetArgvEx), Cvoid,
                                         (Cint, Ref{Ptr{Cvoid}}, Cint),
                                         1, pointer_from_objref(ref0), 0)
        else
            ccall(@pysym(:PySys_SetArgvEx), Cvoid, (Cint, Ptr{Cwstring}, Cint),
                  1, [""], 0)
        end

        # Some Python code checks sys.ps1 to see if it is running
        # interactively, and refuses to be interactive otherwise.
        # (e.g. Matplotlib: see PyPlot#79)
        if isinteractive()
            let sys = pyimport("sys")
                if !haskey(sys, "ps1")
                    sys["ps1"] = ">>> "
                end
            end
        end
    end
end


text:
using IEEE802_15_4

sink       = PacketSink( BPSK, 5, diff_enc = true, verbosity = 0 )
chipfile   = joinpath( dirname(@__FILE__), "ieee802_15_4_900_MHz_sample_chips.dat" )
chipvector = open( readbytes, chipfile )

@time exec( sink, chipvector )


text:
using JuMP
using FlexJuMP
using Gurobi
using LinearAlgebra

# Set the dimensions
n_lines = 2
n_dems = 3

# Setup the uncertainty set parameters
means = [0.; 60.; 10.]
covar = [80. 0 0; 0 80. 0; 0 0 120.]
box_dev = 2 * sqrt.(diag(covar))

# Specify the network details
line_cap = [35, 50]
gen_cap = 100

# Setup the model
m = FlexibilityModel(solver = GurobiSolver(OutputFlag = 0))

# Define variables
@randomvariable(m, d[i = 1:n_dems], mean = means[i])
@recoursevariable(m, a[1:n_lines])
@recoursevariable(m, g)

# Set the line capacity constraints
@constraint(m, [line = 1:n_lines], -line_cap[line] <= a[line])
@constraint(m, [line = 1:n_lines], a[line] <= line_cap[line])

# Set the generator capacity constraints
@constraint(m, 0.0 <= g)
@constraint(m, g <= gen_cap)

# Set the node balance constraints
@constraint(m, a[1] - d[1] == 0)
@constraint(m, -a[1] -a[2] + g - d[2] == 0)
@constraint(m, a[2] - d[3] == 0)

# Define the uncertainty set
setcovariance(m, covar)
setuncertaintyset(m, :Ellipsoid, only_positive = true)
# setuncertaintyset(m, :Hyperbox, [[box_dev]; [box_dev]], only_positive = true)
# setuncertaintyset(m, :PNorm, 1, only_positive = true)
# setuncertaintyset(m, :PNorm, 2, only_positive = true)
# setuncertaintyset(m, :PNorm, Inf, only_positive = true)

# Compute a center to replace the mean if desired
new_mean = findcenteredmean(m, center = :analytic, update_mean = true, only_positive = true)
updated_mean = getmean(m)

# Ensure the mean is valid
if ismeanfeasible(m)
    # Solve
    status = solve(m, U = 1000) # Change slack variable upper bound

    if status == :Optimal
        # Retrieve optimized data
        flexibility_index = getflexibilityindex(m)
        lines = getvalue(a)
        generator = getvalue(g)
        demands = getvalue(d)
        actives = getactiveconstraints(m)

        # Parse data
        data = getflexibilitydata(m)

        # Get the confidence level
        if data.uncertainty_set.name == :Ellipsoid
            conf_lvl = getconfidencelevel(m)
        end

        # Rank the inequality constraints
        rank_data = rankinequalities(m, max_ranks = 4, U = 1000)
    end

    # Estimate the SF index
    SF = findstochasticflexibility(m, use_vulnerability_model = true, use_flexibility_index = true, only_positive = true)
end


text:
using MIToS.MSA
using MIToS.Pfam
using MIToS.PDB
using MIToS.Information
using Printf
using ROCAnalysis

# Download Pfam Stockholm
# msa = downloadpfam('PF08171')

macro output_timings(t,name)
    quote
        @printf "MIToS,%s,%f\n" $name minimum($t)
    end
end

macro timeit(ex,name)
    quote
        t = Float64[]
        for i in 1:6
            e = 1000*(@elapsed $(esc(ex)));
            if i > 1
                # warm up on first iteration
                push!(t, e)
            end
        end
        @output_timings t $name
    end
end

const aln_mapping = read("../data/PF08171.sth", Stockholm, generatemapping=true, useidcoordinates=true);
const col2res = msacolumn2pdbresidue(aln_mapping, "BUB1_YEAST/291-355", "4BL0", "B", "PF08171","../data/4bl0.xml.gz");
const pdb_residues = read("../data/4BL0.xml", PDBML);
const resdict = @residuesdict pdb_residues model "1" chain "B" group "ATOM" residue All;
const cmap = msacontacts(aln_mapping, resdict, col2res);
const ZMIp, MIp = buslje09(aln_mapping, samples=0, clustering=false, lambda=0.0, alphabet=UngappedAlphabet());

@timeit read("../data/PF08171.sth", Stockholm, Matrix{Residue}) "Read Pfam Stockholm MSA"
@timeit read("../data/PF08171.sth", Stockholm) "Read MSA and annotations"
@timeit read("../data/PF08171.sth", Stockholm, generatemapping=true, useidcoordinates=true) "Read MSA and annotations, generate coordinates"
@timeit percentidentity(aln_mapping) "Percent Identity Matrix"
@timeit msacolumn2pdbresidue(aln_mapping, "BUB1_YEAST/291-355", "4BL0", "B", "PF08171","../data/4bl0.xml.gz") "SIFTS residue level mapping"
@timeit read("../data/4BL0.xml", PDBML) "Read PDBML"
@timeit msacontacts(aln_mapping, resdict, col2res) "Protein Contact Map"
@timeit APC!(mapcolpairfreq!(mutual_information, aln_mapping, Counts{Float64,2,GappedAlphabet}(ContingencyTable(Float64,Val{2},GappedAlphabet())))) "Mutual Information APC"
@timeit AUC(MIp, cmap) "AUC (ROC) for contact prediction, MIp"

max_stars_repo_name: diegozea/mitos-benchmarks

max_issues_repo_name: diegozea/mitos-benchmarks

max_forks_repo_name: diegozea/mitos-benchmarks


text:
@testset "Path Methods" begin
    t = 0.001:6π/100:6π
    f = 2 .* max.(0, -(t .- 3 * π)) ./ (t .- 3 * π) .* (sin.(2π * t)) .+ 10 .* max.(0, t .- 3 * π) ./ (t .- 3 * π) .+ max.(0, t .- 3 * π) ./ (t .- 3 * π) .* sin.(4π * t .+ π)
    layers = stFlux((length(f), 1, 1), 2)
    output = layers(f)
    p = pathLocs(1, 1) # an object that we can use as an index into the output
    @test output[p] ≈ output[1][:, 1, 1]
    p = pathLocs(2, (1, 1))
    @test output[p] ≈ output[2][:, 1, 1, 1]
    p = pathLocs(2, (1, 2))
    @test output[p] ≈ output[2][:, 1, 2, 1]
    p = pathLocs(2, (2, 1))
    @test output[p] ≈ output[2][:, 2, 1, 1]
    p = pathLocs(2, (1, :))
    @test output[p] ≈ output[2][:, 1, :, 1]

    # path setindex
    Sx = ScatteredOut((randn(16, 1, 1), randn(11, 32, 1), randn(7, 27, 32, 1)))
    p = pathLocs(0, 5)
    newVal = randn()
    Sx[p] = newVal
    @test Sx[p][1] ≈ Sx.output[1][5, 1, 1]
    @test Sx[p][1] == newVal
    # setting a whole layer
    p = pathLocs(0, :)
    newVal = randn(16)
    Sx[p] = newVal
    @test Sx[p] ≈ Sx.output[1]
    @test Sx[p] ≈ newVal
    # setting a (later) whole layer
    p = pathLocs(2, :)
    newVal = randn(7, 27, 32, 1)
    Sx[p] = newVal
    @test Sx[p] ≈ Sx.output[3]
    @test Sx[p] ≈ newVal

    # setting just a couple of indices
    p = pathLocs(2, (1, 3))
    newVal = randn(7, 1)
    Sx[p] = newVal
    @test Sx[p] ≈ Sx.output[3][:, 1, 3, 1]
    @test Sx[p] ≈ newVal
end


text:
# Codec Interfaces
# ================

"""
### TRANSCODING PROTOCOL (HAPPY PATH)

| NAME              | QTY   | DESC                                              | RETURN                                | HAS DEFAULT?  |
+-------------------+-------+---------------------------------------------------+---------------------------------------+---------------+
| `minoutsize`      | 0:1   |                                                   | minimum output size of `process`      | x             |
| `expectedsize`    | 0:1   |                                                   | expected size of transcoded data      | x             |
| `initialize`      | 0:1   | initialize the codec. called once and only once   |                                       | x             |
| `startproc`       | 0:1   | start processing with the codec                   |                                       | x             |
| `process`         | 1:1   | process data with the codec. called repeatedly    | `(read_size, written_size, :ok)`      | -             |
|                   |       |                                                   | `:end`                                |               |
| `finalize`        | 0:1   | finalize the codec. called once and only once     |                                       | x             |

### ENUMS

MODE
- :idle
- :read
- :write
- :stop
- :close
- :panic

CODE
- :ok
- :end
- :error

### ERROR MGT

todo

"""

abstract type Codec end

minoutsize(codec::Codec, input::Memory)::Int = max(1, div(input.size, 4))
expectedsize(codec::Codec, input::Memory)::Int = input.size
initialize(codec::Codec) = nothing
finalize(codec::Codec)::Nothing = nothing
startproc(codec::Codec, mode::Symbol, error::Error)::Symbol = :ok
process(codec::Codec, input::Memory, output::Memory, error::Error)::Tuple{Int,Int,Symbol} =
    throw(MethodError(process, (codec, input, output, error)))


text:
using SolverPeeker
using Test
using JuMP
using MathOptInterface
const MOI = MathOptInterface
using GLPK

# Start writing tests here, then divide into seperate files

@testset "main" begin

    @testset "lprep.jl" begin

        """ Test not able to conbvert to lp from non-linear constraint. """
        model0 = Model()
        @variable(model0, x >= 0)
        @variable(model0, y >= 0)
        @objective(model0, Max, 3x+5y+2)
        @constraint(model0, x^2+2y*x-y <= 180)

        @test_throws MOI.UnsupportedConstraint get_lpmodel(model0)

        """ Test getting correct A, b and c matrix up to set equality. """
        model1 = Model()
        @variable(model1, x >= 5, Int)
        @variable(model1, y <= 3)
        @variable(model1, 0 <= z <= 50)
        @objective(model1, Max, 3x+5y-2z+2)
        @constraint(model1, x+2y <= 45)
        @constraint(model1, 3z - y >= 4)
        @constraint(model1, 5x+y+2z <= 55)

        lp = get_lpmodel(model1)
        A, b = get_constraint_matrices(lp)
        c = get_objective_vector(lp)
        # From model1, the A matrix, and b vector, when constraints are made less than,
        # sould be equal, up to reshuffeling of elements, to the following:
        mat = [1 2 0; 0 1 -3; 5 1 2]
        vec = [45, -4, 55]
        # Objective only gives variable indices
        obj = [3 5 -2]

        @test issetequal(A, mat)
        @test issetequal(b, vec)
        @test issetequal(c, obj)
        @test lp.objective.constant == 2

        # Add a constraint:
        @constraint(model1, 9x-7y+4-x == 3z)
        lp = get_lpmodel(model1)
        A, b = get_constraint_matrices(lp)
        # The constraint will be split in two, one less than or equal, one greater than or equal,
        # and then the latter will be converted to a less then.
        #All this after summing and moving all the variables to 1 side.
        mat = [mat; 8 -7 -3; -8 7 3]
        vec = [b; -4; 4]

        @test issetequal(A, mat)
        @test issetequal(b, vec)
        # objective should remain the same.
        @test issetequal(c, obj)
        @test lp.objective.constant == 2

        """ Test getting of variables, with correct variable indices. """
        model2 = Model()
        @variable(model2, x >= 5, Int)
        @variable(model2, y <= 3)
        @variable(model2, 0 <= z <= 50)
        @variable(model2, a[1:2], Bin)
        @objective(model2, Max, 3x+5y-2z+sum(a))

        lp2 = get_lpmodel(model2)

        # How the arrays should look like up to set equality
        #greater_than = [5, NaN, 0, NaN, NaN]
        #less_than = [NaN, 3, 50, NaN, NaN]
        #integer = [1, NaN, NaN, NaN, NaN]
        #zero_one = [NaN, NaN, NaN, 1, 1]

        greater_than = get_greater_than_variables(lp2)
        less_than = get_less_than_variables(lp2)
        integer = get_integer_variables(lp2)
        zero_one = get_zero_one_variables(lp2)

        for (index, name) in lp2.var_to_name
            if (name == "x")
                @test greater_than[index.value] == 5
                @test isnan(less_than[index.value])
                @test integer[index.value] == 1
                @test isnan(zero_one[index.value])
            elseif (name == "y")
                @test isnan(greater_than[index.value])
                @test less_than[index.value] == 3
                @test isnan(integer[index.value])
                @test isnan(zero_one[index.value])
            elseif (name == "z")
                @test greater_than[index.value] == 0
                @test less_than[index.value] == 50
                @test isnan(integer[index.value])
                @test isnan(zero_one[index.value])
            elseif (name == "a[1]")
                @test isnan(greater_than[index.value])
                @test isnan(less_than[index.value])
                @test isnan(integer[index.value])
                @test zero_one[index.value] == 1
            elseif (name == "a[2]")
                @test isnan(greater_than[index.value])
                @test isnan(less_than[index.value])
                @test isnan(integer[index.value])
                @test zero_one[index.value] == 1
            else
                # No other names should exist in the model
                @test false
            end
        end

    end

    @testset "recognize.jl" begin

        """ Test difference constraints """

        # Test that only difference constraints return true
        model1 = Model()
        @variable(model1, x>=0, Int)
        @variable(model1, y, Int)
        @variable(model1, z, Int)
        @objective(model1, Min, x-2y-z)
        @constraint(model1, x-y <= 3)
        @constraint(model1, -z+y <= -5)
        @constraint(model1, x-y >= 7)
        @constraint(model1, z-x <= 2)
        @constraint(model1, x-z <= -3)

        lp1 = get_lpmodel(model1)
        @test recognize(lp1, DifferenceConstraints()) == true

        # Add a non-difference constraint, and get false
        @constraint(model1, x+z <= -3)

        lp1 = get_lpmodel(model1)
        @test recognize(lp1, DifferenceConstraints()) == false

        # Test that model with only one variable will return false
        model2 = Model()
        @variable(model2, x)
        @objective(model2, Max, x+2)
        @constraint(model2, x-4 <= 9)

        lp2 = get_lpmodel(model2)
        @test recognize(lp2, DifferenceConstraints()) == false

        # Test if no constraints return false
        model3 = Model()
        @variable(model3, x, Bin)
        @variable(model3, y, Int)
        @variable(model3, z, Int)
        @objective(model3, Max, x+y+z)

        lp3 = get_lpmodel(model3)
        @test recognize(lp3, DifferenceConstraints()) == false

        """ Test AllIntegerVariables """

        @test recognize(lp1, AllIntegerVariables()) == true
        @test recognize(lp2, AllIntegerVariables()) == false
        @test recognize(lp3, AllIntegerVariables()) == false

        """ Test AllIntegerConstraintBounds """

        @test recognize(lp1, AllIntegerConstraintBounds()) == true
        @test recognize(lp2, AllIntegerConstraintBounds()) == true

        model4 = Model()
        @variable(model4, x)
        @variable(model4, y, Int)
        @variable(model4, z, Int)
        @objective(model4, Max, x+y+z)
        @constraint(model4, x-z <= 4.6)

        lp4 = get_lpmodel(model4)
        @test recognize(lp4, AllIntegerConstraintBounds()) == false

    end

    @testset "solve.jl" begin

        model1 = Model()
        @variable(model1, x1, Int)
        @variable(model1, x2, Int)
        @variable(model1, x3, Int)
        @variable(model1, x4, Int)
        @variable(model1, x5, Int)
        @objective(model1, Min, x1+x2-5x3+2x4+6x5)
        @constraint(model1, x1-x2 <= 0)
        @constraint(model1, x1-x5 <= -1)
        @constraint(model1, x2-x5 <= 1)
        @constraint(model1, x3-x1 <= 5)
        @constraint(model1, x4-x1 <= 4)
        @constraint(model1, x4-x3 <= -1)
        @constraint(model1, x5-x3 <= -3)
        @constraint(model1, x5-x4 <= -3)

        lp1 = get_lpmodel(model1)

        # TODO: Test return status
        solve!(lp1, model1, ShortestPath())

    end

    @testset "optimize.jl" begin
        # TODO: Attatch status return and test correct
        # Test feasible unbounded
        model1 = Model()
        @variable(model1, x1, Int)
        @variable(model1, x2, Int)
        @variable(model1, x3, Int)
        @variable(model1, x4, Int)
        @variable(model1, x5, Int)
        @objective(model1, Min, x1+x2-5x3+2x4+6x5)
        @constraint(model1, x1-x2 <= 0)
        @constraint(model1, x1-x5 <= -1)
        @constraint(model1, x2-x5 <= 1)
        @constraint(model1, x3-x1 <= 5)
        @constraint(model1, x4-x1 <= 4)
        @constraint(model1, x4-x3 <= -1)
        @constraint(model1, x5-x3 <= -3)
        @constraint(model1, x5-x4 <= -3)

        @test SolverPeeker.optimize!(model1) == true

        # Test infeasible
        model2 = Model()
        @variable(model2, x>=0, Int)
        @variable(model2, y, Int)
        @variable(model2, z, Int)
        @objective(model2, Min, x-2y-z)
        @constraint(model2, x-y <= 3)
        @constraint(model2, -z+y <= -5)
        @constraint(model2, x-y >= 7)
        @constraint(model2, z-x <= 2)
        @constraint(model2, x-z <= -3)

        @test SolverPeeker.optimize!(model2) == true

        # TODO: Test not recognizeable, no solver attatched
        model3 = Model()
        @variable(model3, x>=0)
        @variable(model3, y)
        @variable(model3, z)
        @objective(model3, Min, x-2y-z)
        @constraint(model3, x-y <= 3)
        @constraint(model3, -z+y <= -5)
        @constraint(model3, x-y >= 7)
        @constraint(model3, z-x <= 2)
        @constraint(model3, 2x-z <= -3)

        @test SolverPeeker.optimize!(model3) == false

        # TODO: Test optimizer is used
        set_optimizer(model3, GLPK.Optimizer)

        @test SolverPeeker.optimize!(model3) == false

    end

end

#include("lprep.jl")
#include("recognize.jl")

max_stars_repo_name: henriean/solver-peeker

max_issues_repo_name: henriean/solver-peeker

max_forks_repo_name: henriean/solver-peeker


text:
# # Defining your own types
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__notebooks/05_DefiningTypes.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__notebooks/05_DefiningTypes.ipynb)
#
# Julia user-defined types are as fast as built-in types. They are easy to
# define:
#

abstract type AbstractFeature end

#-

struct Exon{T} <: AbstractFeature where {T <: Integer}
	start::T
	stop::T
end

#-

Exon(5, 10)

#-

mutable struct MessengerRNA
	sequence::Vector{Char}
	exons::Vector{Exon}
	present::Vector{Bool}
end

#-

mRNA = MessengerRNA(collect("ACTGTTGCATTGCAATTTAAGCAATGGCAAATAACATA"),
					[Exon(5,10), Exon(20,28), Exon(30,34)],
					[true, false, true])

# You can access type fields using a dot:

mRNA.exons

# #### Exercise 1
#
# Define a function that takes a MessengerRNA object and return a list with the
# sequences (as strings) of the exons present in the transcript, i.e.
# `["TTGCAT", "AATAA"]`

## ... your function here ...

max_stars_repo_name: diegozea/JuliaForBioinformatics

max_issues_repo_name: diegozea/JuliaForBioinformatics

max_forks_repo_name: diegozea/JuliaForBioinformatics


text:
# what I figured out for performance.
# send type information T to getindexvalue.
# use genearted @loops instead of arr[index...]

# each allocating a large size array.
# in create_by_array_inner
#  @time subarrays_buf = Array(NTuple{$N,Int}, length(indices))
# in getcondition
#  @time indices::Array{NTuple{$N,Int},1} = Array(NTuple{$N,Int}, length(t))
# in create_bymap_inner
#  @time indexvec::Vector{Int} = Array(Int, lenbyvec)

using Base.Cartesian

# used to test if @inbounds results in any substantial speed up.
#macro debuginbounds(x)
#  esc(x)
#  #:(@inbounds($(esc(x))))
#end

"""

`selectfunc(t::LabeledArray, c, b, a)`

main select function. This function is internal and is meant to be used via `selct`.
`selectfunc` takes a table `t`, condition `c`, aggreagtion by rule `b`, aggregate function `a`.

* t : a LabeledArray
* c : an array of functions `(t, inds)` -> boolean array. `t` is a LabeledArray and inds is either `nothing` to choose the entire `t` or an array of Cartesian indices.
* b : an array of arrays of pairs of field name => function `(t, inds)` -> array. `t` is a LabeledArray and inds is either `nothing` to choose the entire `t` or an array of Cartesian indices.
* c : an array of pairs of field name => function `(t, inds)` -> array. `t` is a LabeledArray and inds is either `nothing` to choose the entire `t` or an array of Cartesian indices.

"""
selectfunc{N}(t::LabeledArray{TypeVar(:T),N}, c, b, a) = begin
  nocond = isempty(c)
  noby = isempty(b)
  noagg = isempty(a)
  if nocond
    if noby
      if noagg
        t
      else
        if isa(t.data, DictArray)
          resa = [(k, onea(t,nothing)) for (k,onea) in a]
          if all(x->!isa(x[2], AbstractArray), resa)
            LDict(resa...)
          else
            LabeledArray(darr(resa...), t.axes)
          end
        else
          error("this should not happen: the data part of the LabeledArray is not a DictArray, but agg expression is provided.")
        end
      end
    else
      # nocond, yesby
      selected_cartesian_indices_all::Vector{NTuple{N,Int}} = getcondition(t, [(tbl, inds) -> inds])
      byvecs = [darr([(k, onebcomp(t, selected_cartesian_indices_all)) for (k,onebcomp) in oneb]...) for oneb in b]
      bymaps = [create_bymap(byvec) for byvec in byvecs]
      create_by_array(t, a, selected_cartesian_indices_all, byvecs, bymaps)
    end
  else
    # yescond
    selected_cartesian_indices::Vector{NTuple{N,Int}} = getcondition(t, c)
    if noby
      # results of a are not vectors.
      if !noagg && isa(t.data, DictArray)
        resa = [(k, onea(t,selected_cartesian_indices)) for (k,onea) in a]
        if all(x->!isa(x[2], AbstractArray), resa)
          return LDict(resa...)
        end
      end

      # yescond, noby
      aggvecs = if noagg
        if isa(t.data, DictArray)
          darr([(k, select_cartesian_elements_helper(onefld, selected_cartesian_indices)) for (k,onefld) in t.data.data]...)
          #above gives better result than below.
          #darr([(k, eltype(onefld)[onefld[i...] for i in selected_cartesian_indices]) for (k,onefld) in t.data.data]...)
        else
          tdata = t.data
          select_cartesian_elements_helper(tdata, selected_cartesian_indices)
          #above gives better result than below.
          #eltype(tdata)[tdata[i...] for i in selected_cartesian_indices]
        end
      else
        if isa(t.data, DictArray)
          darr([(k, onea(t, selected_cartesian_indices)) for (k,onea) in a]...)
        else
          error("this should not happen: the data part of the LabeledArray is not a DictArray, but agg expression is provided.")
        end
      end
      create_square_array(aggvecs, t, selected_cartesian_indices)
    else
      # yescond, yesby
      byvecs = [darr([(k, onebcomp(t, selected_cartesian_indices)) for (k,onebcomp) in oneb]...) for oneb in b]
      bymaps = [create_bymap(byvec) for byvec in byvecs]
      create_by_array(t, a, selected_cartesian_indices, byvecs, bymaps)
    end
  end
end

#select_cartesian_elements_helper{T,N}(fld::AbstractArray{T,N}, indices::Vector{NTuple{N,Int}}) = begin
#  result = Array(T, length(indices))
#  select_cartesian_elements_helper2(fld

@generated select_cartesian_elements_helper{T,N}(fld::AbstractArray{T,N}, indices::Vector{NTuple{N,Int}}) = quote
  result = Array(T, length(indices))
  #result::T = similar(fld, length(indices))
  for i in eachindex(result)
    @inbounds index = indices[i]
    @inbounds result[i] = @nref $N fld d->index[d]
  end
  result
end

@generated select_cartesian_elements_helper{T<:AbstractFloat,N,A}(fld::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}},
                                                                indices::Vector{NTuple{N,Int}}) = quote
  result = Array(T, length(indices))
  #result::T = similar(fld, length(indices))
  fldunderlying = fld.a.data
  for i in eachindex(result)
    @inbounds index = indices[i]
    @inbounds result[i] = @nref $N fldunderlying d->index[d]
  end
  AbstractArrayWrapper(FloatNAArray(result))
end

"""

`updatefunc(t::LabeledArray, c, b, a)`

main update function. This function is internal and is meant to be used via `update`.
`updatefunc` takes a table `t`, condition `c`, aggreagtion by rule `b`, aggregate function `a`.

* `t` : a LabeledArray
* `c` : an array of functions `(t, inds)` -> boolean array. `t` is a LabeledArray and inds is either `nothing` to choose the entire `t` or an array of Cartesian indices.
* `b` : an array of arrays of pairs of field name => function `(t, inds)` -> array. `t` is a LabeledArray and inds is either `nothing` to choose the entire `t` or an array of Cartesian indices.
* `c` : an array of pairs of field name => function `(t, inds)` -> array. `t` is a LabeledArray and inds is either `nothing` to choose the entire `t` or an array of Cartesian indices.

"""
updatefunc{N}(t::LabeledArray{TypeVar(:T),N}, c, b, a) = begin
  nocond = isempty(c)
  noby = isempty(b)
  noagg = isempty(a)
  if noagg
    t
  elseif nocond
    # agg but nocond noby
    if noby
      addres = map(a) do konea
        konea[1] => konea[2](t, nothing)
      end
      LabeledArray(darr(merge(t.data.data, LDict(addres...))...), t.axes)
    else
      # nocond, yesby
      selected_cartesian_indices_all::Vector{NTuple{N,Int}} = getcondition(t, [(tbl, inds) -> inds])
      flattened_b = [b...;]
      byvec = darr([(k, onebcomp(t, selected_cartesian_indices_all)) for (k,onebcomp) in flattened_b]...)
      bymap = create_bymap(byvec, true)
      update_by_array(t, a, selected_cartesian_indices_all, byvec, bymap)
    end
  else
    # agg with cond
    selected_cartesian_indices::AbstractArray{NTuple{N,Int},1} = getcondition(t, c)
    if isempty(selected_cartesian_indices)
      return t
    end
    if noby
      result = copy(t.data.data)
      addres = map(a) do konea
        k = konea[1]
        onea = konea[2]
        aggvec = onea(t, selected_cartesian_indices)
        if haskey(result, k)
          resultk = copy(result[k])
          i = 1
          for coords in selected_cartesian_indices
            resultk[coords...] = aggvec[i]
            i += 1
          end
        else
          resultk = similar(aggvec, size(t.data))
          setna!(resultk, :)
          i = 1
          for coords in selected_cartesian_indices
            resultk[coords...] = aggvec[i]
            i += 1
          end
        end
        k => resultk
      end
      LabeledArray(darr(merge(result, LDict(addres...))...), t.axes)
    else
      # yes agg, by, cond
      flattened_b = [b...;]
      byvec = darr([(k, onebcomp(t, selected_cartesian_indices)) for (k,onebcomp) in flattened_b]...)
      bymap = create_bymap(byvec, true)
      update_by_array(t, a, selected_cartesian_indices, byvec, bymap)
    end
  end
end

# creates a by array for selectfunc.
create_by_array(t, a, indices, byvecs, bymaps) = begin
  create_by_array_inner(t, a, indices, byvecs, bymaps, Array(Int,ntuple(_->0, length(bymaps))...))
end

# N is the original array dimensions.
# M is the byarray dimentions.
@generated create_by_array_inner{N,M}(t::AbstractArray{TypeVar(:T),N}, a, indices::Vector{NTuple{N,Int}}, byvecs, bymaps, dummy::Array{Int,M}) = quote
  labels::NTuple{$M} = (map(x->x[1], bymaps)...)
  # assumes there is at least one label for each axis.
  ressize::NTuple{$M,Int} = (map(x->length(x[1][1]), bymaps)...)
  indexvecs::NTuple{$M,Vector{Int}} = (map(x->x[2], bymaps)...)
  indcountmat::Array{Int,$M} = fill(0, ressize)
  for i in eachindex(indices)
    @inbounds @nexprs $M j->(ids_j = indexvecs[j][i])
    @inbounds @nref($M,indcountmat,ids) += 1
  end
  subarrays_buf = Array(NTuple{$N,Int}, length(indices))
  index_subarrays = similar(indcountmat, SubArray{NTuple{$N,Int},1,Array{NTuple{$N,Int},1},Tuple{UnitRange{Int}},1})
  indsumsofar::Int = 0
  lenindcountmat::Int = length(indcountmat)
  for i in 1:lenindcountmat
    prevsum = indsumsofar
    @inbounds indsumsofar += indcountmat[i]
    @inbounds index_subarrays[i] = sub(subarrays_buf, prevsum+1:indsumsofar)
  end
  # change the meaning of indcountmat, not to create another buffer.
  fill!(indcountmat, 0)
  counter::Int = 1
  for i in eachindex(indices)
    @inbounds @nexprs $M j->(ids_j = indexvecs[j][i])
    @inbounds current_count = @nref($M,indcountmat,ids) + 1
    @inbounds @nref($M,indcountmat,ids) = current_count
    @inbounds subarray = @nref($M, index_subarrays, ids)
    @inbounds subarray[current_count] = indices[i]
  end
  effa = if isempty(a)
    [(k, ((t,r) -> selectfield(t,k,r))) for (k,_) in t.data.data]
  else
    a
  end

  # There is a problem when an element in index_subarrays is empty. In that case, the type is not properly marked.
  # This is ultimately due to strange map function behavior. For example,
  #   map((x,y)->x==y, Int[], Int[])
  # is an Int array, not a boolean array (Granted that I couldn't see a better way, though).
  # So, comment out this:
  # DictArray(LDict([(k, map(r->lift_naelem(onea(t,r)), index_subarrays)) for (k, onea) in a]...))
  # and implement another method, by creating an empty array for each onea in a, and fill it up.
  # But first, we need to check the type of each element in each onea, which can be done by examining any
  # non zero length array results. We assume that all elements have the same type.
  restypes = map(effa) do konea
    onea = konea[2]
    for inds in index_subarrays
      # perhaps it is okay to just use one index.
      # but do not attempt it now.
      #return typeof(apply_nullable(onea(t,inds[1:1])))
      if !isempty(inds)
        return typeof(apply_nullable(onea(t,inds)))
      end
    end
    return Any
  end
  # Now, let's create empty arrays to fill up results later.
  results = map(tp -> similar(index_subarrays, tp), restypes)
  for (konea, result) in zip(effa, results)
    fill_byarrays!(result, konea[2], index_subarrays, t)
  end
  aggvecs = DictArray(LDict([(konea[1], result) for (konea, result) in zip(effa, results)]...))
  newaxes = map(byvecs, labels) do byvec, label
    DictArray(byvec.data.keys, collect(label))
  end
  LabeledArray(aggvecs, (newaxes...))
end

@generated fill_byarrays!{T,M}(result::AbstractArray{Nullable{T},M}, onea, index_subarrays::AbstractArray{TypeVar(:V),M}, t::LabeledArray) = quote
  @nloops $M i result begin
    inds = @nref($M,index_subarrays,i)
    value = if isempty(inds)
      Nullable{T}()
    else
      apply_nullable(onea(t, inds))
    end
    @nref($M,result,i) = value
  end
end

bymap_coords_calculate_helper{T}(bymap::Dict{T,Int}, byvec, bytype::Type{T}, i::Int) = begin
  v::T = getindexvalue(byvec,bytype,i)
  bymap[v]::Int
end

create_bymap(byvec::DictArray, skip_ordering::Bool=false) = create_bymap_inner(byvec, Tuple{[eltype(v) for v in byvec.data.values]...}, skip_ordering)

create_bymap_inner{T}(byvec::DictArray, ::Type{T}, skip_ordering::Bool) = if isa(byvec, DictArray)
  # resmap: (d dim tuple) for the key tuple => the coordiante in the result byarray.
  lenbyvec::Int = length(byvec)
  resmap = Dict{T,Int}()
  # N 1 dimensional vectors to store the result axis labels.
  labelvecs = ntuple(length(byvec.data)) do d
    eltype(byvec.data.values[d])[]
  end
  # well, this type annotation is necessary to enhance performance!
  # don't know why the type is not inferred stably.
  indexvec::Vector{Int} = Array(Int, lenbyvec)
  counter = 1
  for i in 1:lenbyvec
    value::T = getindexvalue(byvec, T, i)
    # index is the index in the output byarray along some direction.
    # why is this more efficient than get! ?
    index::Int = if haskey(resmap, value)
      resmap[value]
    else
      push_labelvecs_inner!(labelvecs, value)
      tempcounter = counter
      counter += 1
      resmap[value] = tempcounter
      tempcounter
    end
    @inbounds indexvec[i] = index
  end
  # check to see if reordeingr indices is necessary.
  tosort_indices = Bool[method_exists(isless,(peel_nullabletype(eltype(x)),peel_nullabletype(eltype(x)))) for x in labelvecs]
  nochange_flag = skip_ordering || !any(tosort_indices)
  if nochange_flag
    labelvecs, indexvec
  else
    # reorder indices.
    toorder_labels = collect(zip(labelvecs[tosort_indices]...))
    neworder = sortperm(toorder_labels, lt=byarray_isless)
    len = length(toorder_labels)
    # well, this type annotation is necessary to enhance performance!
    # don't know why the type is not inferred correctly.
    inverse_order::Vector{Int} = Array(Int, len)
    for i in eachindex(neworder)
      @inbounds inverse_order[neworder[i]] = i
    end
    for labelvec in labelvecs
      permute!(labelvec, neworder)
    end
    for i in eachindex(indexvec)
      @inbounds indexvec[i] = inverse_order[indexvec[i]]
    end
    labelvecs, indexvec
  end
else
  error("did not implelment yet.")
end

# this does not gain any so far.
push_labelvecs_inner!{T1}(labelvecs::Tuple{Vector{T1}}, value::Tuple{T1}) = begin
  push!(labelvecs[1], value[1])
end
push_labelvecs_inner!{T1,T2}(labelvecs::Tuple{Vector{T1},Vector{T2}}, value::Tuple{T1,T2}) = begin
  push!(labelvecs[1], value[1])
  push!(labelvecs[2], value[2])
end
push_labelvecs_inner!{T1,T2,T3}(labelvecs::Tuple{Vector{T1},
                                                 Vector{T2},
                                                 Vector{T3}},
                                value::Tuple{T1,T2,T3}) = begin
  push!(labelvecs[1], value[1])
  push!(labelvecs[2], value[2])
  push!(labelvecs[3], value[3])
end
push_labelvecs_inner!{T1,T2,T3,T4}(labelvecs::Tuple{Vector{T1},
                                                   Vector{T2},
                                                   Vector{T3},
                                                   Vector{T4}},
                                   value::Tuple{T1,T2,T3,T4}) = begin
  push!(labelvecs[1], value[1])
  push!(labelvecs[2], value[2])
  push!(labelvecs[3], value[3])
  push!(labelvecs[4], value[4])
end
push_labelvecs_inner!(labelvecs::Tuple, value::Tuple) = begin
  for (v,labelvec) in zip(value, labelvecs)
    push!(labelvec, v)
  end
end

peel_nullabletype{T}(::Type{Nullable{T}}) = T
peel_nullabletype(::Type{Nullable}) = Any

byarray_isless{T}(x::Nullable{T}, y::Nullable{T}) =
  (x.isnull && !y.isnull) || (!x.isnull && !y.isnull && isless(x.value, y.value))
byarray_isless{T1,T2}(x::Tuple{Nullable{T1},
                               Nullable{T2}},
                      y::Tuple{Nullable{T1},
                               Nullable{T2}}) = begin
  if byarray_isless(x[1], y[1])
    true
  elseif byarray_isless(y[1], x[1])
    false
  else
    byarray_isless(x[2], y[2])
  end
end
byarray_isless{T1,T2,T3}(x::Tuple{Nullable{T1},
                                  Nullable{T2},
                                  Nullable{T3}},
                         y::Tuple{Nullable{T1},
                                  Nullable{T2},
                                  Nullable{T3}}) = begin
  if byarray_isless(x[1], y[1])
    true
  elseif byarray_isless(y[1], x[1])
    false
  elseif byarray_isless(x[2], y[2])
    true
  elseif byarray_isless(y[2], x[2])
    false
  else
    byarray_isless(x[3], y[3])
  end
end
byarray_isless{T1,T2,T3,T4}(x::Tuple{Nullable{T1},
                                  Nullable{T2},
                                  Nullable{T3},
                                  Nullable{T4}},
                         y::Tuple{Nullable{T1},
                                  Nullable{T2},
                                  Nullable{T3},
                                  Nullable{T4}}) = begin
  if byarray_isless(x[1], y[1])
    true
  elseif byarray_isless(y[1], x[1])
    false
  elseif byarray_isless(x[2], y[2])
    true
  elseif byarray_isless(y[2], x[2])
    false
  elseif byarray_isless(x[3], y[3])
    true
  elseif byarray_isless(y[3], x[3])
    false
  else
    byarray_isless(x[4], y[4])
  end
end
byarray_isless{N}(x::NTuple{N}, y::NTuple{N}) = begin
  for i in 1:N
    if byarray_isless(x[i], y[i])
      return true
    elseif byarray_isless(y[i], x[i])
      return false
    end
  end
  false
end

create_square_array{T,N}(aggvec, t::AbstractArray{T,N}, inds::Vector{NTuple{N,Int}}) = begin
  sz::NTuple{N,Int} = size(t)
  boolvecs::NTuple{N,BitVector} = ntuple(d->falses(sz[d]), N)
  create_square_array_helper1!(boolvecs, inds)
  wherevecs::NTuple{N,Vector{Int}} = ntuple(d->find(boolvecs[d]), N)
  reversemapvec::NTuple{N,Vector{Int}} = map(boolvecs, wherevecs) do boolvec, wherevec
    v = Array(Int, length(boolvec))
    @inbounds v[wherevec] = 1:length(wherevec)
    v
  end
  ressize = map(length, wherevecs)
  result = similar(aggvec, ressize)

  setna!(result)
  newinds = create_square_array_helper2(inds, reversemapvec)
  setcartesian!(result, aggvec, newinds)
  # result is the data field of the return value. Now work on the axes part.
  @inbounds newaxes = map((axis, wherevec) -> isa(axis, DefaultAxis) ? DefaultAxis(length(wherevec)) : axis[wherevec], t.axes, wherevecs)
  LabeledArray(result, (newaxes...))
end

@generated create_square_array_helper1!{N}(boolvecs::NTuple{N,BitVector}, inds::Vector{NTuple{N,Int}}) = quote
  for ind in inds
    @nexprs $N d->boolvecs[d][ind[d]] = true
  end
end

@generated create_square_array_helper2{N}(inds::AbstractVector{NTuple{N,Int}}, reversemapvec::NTuple{N,Vector{Int}}) = quote
  newinds = Array(NTuple{$N,Int}, length(inds))
  for i in eachindex(newinds)
    tpl = @ntuple $N d->reversemapvec[d][inds[i][d]]
    @inbounds newinds[i] = tpl
  end
  newinds
end
# the current version of create_square_array_helper in julia 0.4.1 turns out to be definitely faster than below.
# map(ind->ntuple(d->reversemapvec[d][ind[d]], N), inds)

setcartesian!{N}(tgt::DictArray, src::DictArray, inds::AbstractArray{NTuple{N,Int},1}) = begin
  for (k, v) in src.data
    @inbounds onefield = tgt.data[k]
    setcartesian_inner!(inds, onefield, v)
  end
  tgt
end
setcartesian!{T,N}(tgt::AbstractArray{Nullable{T},N}, src::AbstractArray{Nullable{T},1}, inds::AbstractArray{NTuple{N,Int},1}) = begin
  setcartesian_inner!(inds, tgt, src)
  tgt
end

setcartesian_inner!{T,N}(inds::AbstractArray{NTuple{N,Int},1}, onefield::AbstractArray{T,N}, v::AbstractArray{T,1}) = begin
  for (ind,onev) in zip(inds,v)
    @inbounds onefield[ind...] = onev
  end
  # this hardly gains any.
  #for i in eachindex(inds)
  #  @inbounds onefield[inds[i]...] = v[i]
  #end
end

# A view of an array, to be used in selectfield below.
# It won't be too difficult to implement because we assume it is a one dimensional array,
# determined by a vector of cartesian indices.
immutable SubArrayView{T,N,A,I} <: AbstractVector{T}
  data::A
  indices::I
end

# a syb array view of an array to select a portion of a `LabeledArray` during `selct` or `update`.
SubArrayView{T,N}(data::AbstractArray{T,N}, indices::AbstractVector{NTuple{N,Int}}) = SubArrayView{T,N,typeof(data),typeof(indices)}(data, indices)
SubArrayView{T,N,A}(data::AbstractArrayWrapper{T,N,A}, indices::AbstractVector{NTuple{N,Int}}) = SubArrayView{T,N,A,typeof(indices)}(data.a, indices)
SubArrayView{T,N,A}(data::AbstractArrayWrapper{Nullable{T},N,FloatNAArray{T,N,A}}, indices::AbstractVector{NTuple{N,Int}}) = FloatNAArray(SubArrayView{T,N,A,typeof(indices)}(data.a.data, indices))
Base.getindex{T,N,A,I}(arr::SubArrayView{T,N,A,I}, arg::Int) = arr.data[arr.indices[arg]...]
# Let's make the array immutable at this stage.
#Base.setindex!{T,N}(arr::SubArrayView{T,N}, v::T, arg::Int) = setindex!(arr.data, v, arr.indices[arg]...)
Base.start(iter::SubArrayView) = 1
Base.next{T,N,A,I}(iter::SubArrayView{T,N,A,I}, state) = (iter.data[iter.indices[state]...], state+1)
Base.done(iter::SubArrayView, state) = state > length(iter.indices)
Base.size(arr::SubArrayView) = (length(arr.indices),)
Base.eltype{T,N,A,I}(::Type{SubArrayView{T,N,A,I}}) = T
Base.endof(arr::SubArrayView) = length(arr.indices)
Base.linearindexing{T,N,A,I}(::Type{SubArrayView{T,N,A,I}}) = Base.LinearFast()

"""

`selectfield(t, fld, inds)` : select a field whose name is `fld` at cartesian coordinates `inds` in a LabeledArray `t`.
If `inds` is `nothing`, it chooses an entire `fld` from `t`.

"""
function selectfield end

selectfield(t::LabeledArray, fld, ::Void) = begin
  onefld = selectfield(t, fld)
  wrap_array(onefld)
end
selectfield(t::LabeledArray, fld, inds) = AbstractArrayWrapper(SubArrayView(selectfield(t,fld), inds))
selectfield{K,N}(t::LabeledArray, fld::AbstractArray{Nullable{K},N}, ::Void) = begin
  # use a slower version: find a field for each ind.
  @inbounds r=simplify_array([onefld.isnull ? Nullable{Any}() : selectfield(t,onefld.value)[i] for (i,onefld) in zip(eachindex(fld), fld)])
  reshape(r, size(fld))
end

selectfield{K,N}(t::LabeledArray, fld::AbstractVector{Nullable{K}}, inds::AbstractArray{NTuple{N,Int},1}) = begin
  if length(fld) != length(inds)
    throw(ArgumentError("key lengths do not match with indices length."))
  end
  # use a slower version: find a field for each ind.
  @inbounds r=simplify_array([onefld.isnull ? Nullable{Any}() : selectfield(t,onefld.value)[i...] for (i,onefld) in zip(inds,fld)])
  r
end

# create an array of cartesian indices of a LabeledArray that satisfies
# the conditions in c.
@generated getcondition{T,N}(t::LabeledArray{T,N}, c) = quote
  indices::Array{NTuple{$N,Int},1} = Array(NTuple{$N,Int}, length(t))
  index = 1
  @nloops $N i t begin
    tupleind::NTuple{$N,Int} = @ntuple $N i
    @inbounds indices[index] = tupleind
    index += 1
  end
  foldl((inds, onec) -> onec(t, inds), indices, c)
end

# larray -> indices -> field name -> element.
# Implementationwise, it is easier to do larray -> fieldname -> indices -> element.
# This and other *LazilySelectedLabeledArray class help to use the other mapping above.
immutable LazilySelectedLabeledArray{T,N,AXES<:Tuple,TN,I} <: Associative
  larray::LabeledArray{T,N,AXES,TN}
  indices::I #::{Array{NTuple{N,Int}}
end

immutable IGNALazilySelectedLabeledArray{T,N,AXES<:Tuple,TN,I} <: Associative
  larray::LabeledArray{T,N,AXES,TN}
  indices::I #::Array{NTuple{N,Int}}
end

immutable ISNALazilySelectedLabeledArray{T,N,AXES<:Tuple,TN,I} <: Associative
  larray::LabeledArray{T,N,AXES,TN}
  indices::I #Array{NTuple{N,Int}}
end

Base.get{K}(d::LazilySelectedLabeledArray, k::AbstractArray{Nullable{K}}, _) = selectfield(d.larray, k, d.indices)
Base.get{K}(d::IGNALazilySelectedLabeledArray, k::AbstractArray{Nullable{K}}, _) = igna(selectfield(d.larray, k, d.indices))
Base.get{K}(d::ISNALazilySelectedLabeledArray, k::AbstractArray{Nullable{K}}, _) = isna(selectfield(d.larray, k, d.indices))
Base.get{K}(d::LazilySelectedLabeledArray, k::K, _) = selectfield(d.larray, k, d.indices)
Base.get{K}(d::IGNALazilySelectedLabeledArray, k::K, _) = igna(selectfield(d.larray, k, d.indices))
Base.get{K}(d::ISNALazilySelectedLabeledArray, k::K, _) = isna(selectfield(d.larray, k, d.indices))
Base.length{T,N,AXES<:Tuple,TN}(d::LazilySelectedLabeledArray{T,N,AXES,TN,Void}) = length(d.larray)
Base.length{T,N,AXES<:Tuple,TN,I}(d::LazilySelectedLabeledArray{T,N,AXES,TN,I}) = length(d.indices)
Base.length{T,N,AXES<:Tuple,TN}(d::IGNALazilySelectedLabeledArray{T,N,AXES,TN,Void}) = length(d.larray)
Base.length{T,N,AXES<:Tuple,TN,I}(d::IGNALazilySelectedLabeledArray{T,N,AXES,TN,I}) = length(d.indices)
Base.length{T,N,AXES<:Tuple,TN}(d::ISNALazilySelectedLabeledArray{T,N,AXES,TN,Void}) = length(d.larray)
Base.length{T,N,AXES<:Tuple,TN,I}(d::ISNALazilySelectedLabeledArray{T,N,AXES,TN,I}) = length(d.indices)
Base.size{T,N,AXES<:Tuple,TN}(d::LazilySelectedLabeledArray{T,N,AXES,TN,Void}) = size(d.larray)
Base.size{T,N,AXES<:Tuple,TN,I}(d::LazilySelectedLabeledArray{T,N,AXES,TN,I}) = size(d.indices)
Base.size{T,N,AXES<:Tuple,TN}(d::ISNALazilySelectedLabeledArray{T,N,AXES,TN,Void}) = size(d.larray)
Base.size{T,N,AXES<:Tuple,TN,I}(d::ISNALazilySelectedLabeledArray{T,N,AXES,TN,I}) = size(d.indices)
Base.size{T,N,AXES<:Tuple,TN}(d::IGNALazilySelectedLabeledArray{T,N,AXES,TN,Void}) = size(d.larray)
Base.size{T,N,AXES<:Tuple,TN,I}(d::IGNALazilySelectedLabeledArray{T,N,AXES,TN,I}) = size(d.indices)
Base.similar(d::LazilySelectedLabeledArray, element_type::Type, dim, dims...) = similar(d.larray, element_type, dim, dims...)
Base.similar(d::IGNALazilySelectedLabeledArray, element_type::Type, dim, dims...) = similar(d.larray, element_type, dim, dims...)
Base.similar(d::ISNALazilySelectedLabeledArray, element_type::Type, dim, dims...) = similar(d.larray, element_type, dim, dims...)

igna(d::LazilySelectedLabeledArray) = IGNALazilySelectedLabeledArray(d.larray, d.indices)
isna(d::LazilySelectedLabeledArray) = ISNALazilySelectedLabeledArray(d.larray, d.indices)

"""

`sel(func, t [; c=[], b=[], a=[]])`

an intermediate `select`/`update` function to connect `selectfunc`/`updatefunc` and `@select`(`selct`)/`@update`(`update`).

##### Arguments

Below `t'` is an object such that `t'[k]` for a field name `k` gives the corresponding array
for the field `k` in the table `t` at the coordinates selected so far.
* `func` : selectfunc or updatefunc.
* `t` : a LabeledArray.
* `c` : an array of conditions of type `t'` -> nullable boolean array.
* `b` : an array of arrays of pairs from field names to by functions specified as `t'` -> a nullable array.
* `a` : an array of pairs from field names to aggregate functions specified as `t'` -> a nullable array.

"""
sel(func, t; c=[], b=[], a=[]) = begin
  efft = isa(t, DictArray) ? LabeledArray(t) : t
  cfuncs = map(c) do onec
    if isa(onec, Function)
      (d, inds) -> begin
        @inbounds r = inds[ignabool(onec(LazilySelectedLabeledArray(efft, inds)))]
        r
      end
    else
      (d, inds) -> begin
        @inbounds r = inds[ignabool(selectfield(d, onec, inds))]
        r
      end
    end
  end
  bfuncs = map(b) do oneb
    map(oneb) do subb
      if !isa(subb, Pair) && !isa(subb, Tuple)
        subb => (d, inds) -> selectfield(d, subb, inds)
      elseif isa(subb[2], Function)
        subb[1] => (d, inds) -> subb[2](LazilySelectedLabeledArray(efft, inds))
      else
        subb[1] => (d, inds) -> selectfield(d, subb[2], inds)
      end
    end
  end
  afuncs = map(a) do onea
    if !isa(onea, Pair) && !isa(onea, Tuple)
      onea => (d, inds) -> selectfield(d, onea, inds)
    elseif isa(onea[2], Function)
      onea[1] => (d, inds) -> onea[2](LazilySelectedLabeledArray(efft, inds))
    else
      onea[1] => (d, inds) -> selectfield(d, onea[2], inds)
    end
  end
  if isa(t, DictArray)
    result = func(efft, cfuncs, bfuncs, afuncs)
    if isa(result, LabeledArray) && all(map(axis->isa(axis,DefaultAxis), result.axes))
      result.data
    else
      result
    end
  else
    func(efft, cfuncs, bfuncs, afuncs)
  end
end

"""

`@select(t, args...)`

Select macro transforms a `LabeledArray` or `DictArray` into another by choosing / grouping / aggregating.

##### Arguments

Below `t'` is an object such that `t'[k]` for a field name `k` gives the corresponding array
for the field `k` in the table `t` at the coordinates selected so far.

* `t` : a `LabeledArray` or `DictArray`.
* `args...` : each argument can be field names, fieldname=>(t'->nullable array function) pair, `by[...]` or `where[...]`. A fieldname=function can be used instead of the pair notation if the fieldname is a symbol. The first `by[...]` has an array of similar expressions and determines the 1st axis. The second `by[...]` similarly determines the 2nd axis. The output `LabeledArray` will have dimensions of the number of `by[...]` clauses, or the original dimensions if no `by[...]` is provided. `where[...]` has (`t'`->nullable boolean array) functions inside and chooses the appropriate portion of the original array. Multiple `where[...]` will simply be combined.

##### Function Specification

Note that a function (`t'`->nullable array) is expressed by some **expression** with variable names with underscores. The **expression** is converted into `t''->**expression**. Symbols with underscores are converted in the following way:

* `_k` : translates to `t'[k]`. The field name `k` should be a symbol.
* `_!k` : translates to `isna(t')[k]`. It gives a boolean array to denote whether an element is null. The field name `k` should be a symbol.
* `_` : translates to `t'` if none of the above is applicable. It is useful when the field name is not a symbol.

##### Return

A `LabeledArray` transformed by `args...` if `t` is a `LabeledArray`.
If `t` is `DictArray` and the transformed `LabeledArray` has `DefaultAxis` along each direction, the return value is also a `DictArray`. Otherwise it is a `LabeledArray`.

##### Examples

```julia
julia> t = @larr(a=1:10, b=[1,2,3,NA,NA,NA,1,1,2,3], c=[:x,:x,:x,:x,:y,:y,:y,:z,:z,:z])
10 LabeledArray

   |a  b c
---+-------
1  |1  1 x
2  |2  2 x
3  |3  3 x
4  |4    x
5  |5    y
6  |6    y
7  |7  1 y
8  |8  1 z
9  |9  2 z
10 |10 3 z

julia> @select(t, :a, :b=>_b .* 2)
10 LabeledArray

   |a  b
---+-----
1  |1  2
2  |2  4
3  |3  6
4  |4
5  |5
6  |6
7  |7  2
8  |8  2
9  |9  4
10 |10 6

julia> @select(t, :a, :b=>_b .* 2, where[_c .!= :z], where[_a .> 2])
5 LabeledArray

  |a b
--+----
1 |3 6
2 |4
3 |5
4 |6
5 |7 2

julia> @select(t, a=mean(_a), b=sum(_b .* 2), by[d=_b .* 2], by[:c])
4 x 3 LabeledArray

c |x     |y     |z
--+------+------+-------
d |a   b |a   b |a    b
--+------+------+-------
  |4.0 0 |5.5 0 |
2 |1.0 2 |7.0 2 |8.0  2
4 |2.0 4 |      |9.0  4
6 |3.0 6 |      |10.0 6
```

"""
macro select(t, args...)
  select_inner(:selectfunc, t, args...)
end

"""

`selct(t, agg... [; by=[...]..., where=[...]...])`

Select a `LabeledArray` or `DictArray` into another by choosing / grouping / aggregating.

##### Arguments

Below `t'` is an object such that `t'[k]` for a field name `k` gives the corresponding array
for the field `k` in the table `t` at the coordinates selected so far.

* `t` : a `LabeledArray` or `DictArray`.
* `agg...` : each argument can be field names, fieldname=>(t'->nullable array function) pair, `by[...]` or `where[...]`. A fieldname=function can be used instead of the pair notation if the fieldname is a symbol.
* `by=[...]`: the first `by=[...]` has an array of similar expressions and determines the 1st axis. The second `by=[...]` similarly determines the 2nd axis. The output `LabeledArray` will have dimensions of the number of `by[...]` clauses, or the original dimensions if no `by[...]` is provided.
* `where=[...]`: has (`t'`->nullable boolean array) functions inside and chooses the appropriate portion of the original array. Multiple `where=[...]` will simply be combined.

##### Return

A `LabeledArray` transformed by `args...` if `t` is a `LabeledArray`.
If `t` is `DictArray` and the transformed `LabeledArray` has `DefaultAxis` along each direction, the return value is also a `DictArray`. Otherwise it is a `LabeledArray`.

##### Examples

```julia
julia> t = @larr(a=1:10, b=[1,2,3,NA,NA,NA,1,1,2,3], c=[:x,:x,:x,:x,:y,:y,:y,:z,:z,:z])
10 LabeledArray

   |a  b c
---+-------
1  |1  1 x
2  |2  2 x
3  |3  3 x
4  |4    x
5  |5    y
6  |6    y
7  |7  1 y
8  |8  1 z
9  |9  2 z
10 |10 3 z

julia> selct(t, :a, :b=>d->d[:b] .* 2)
10 LabeledArray

   |a  b
---+-----
1  |1  2
2  |2  4
3  |3  6
4  |4
5  |5
6  |6
7  |7  2
8  |8  2
9  |9  4
10 |10 6

julia> selct(t, :a, :b=>d->d[:b] .* 2, where=[d->d[:c] .!= :z], where=[d->d[:a] .> 2])
5 LabeledArray

  |a b
--+----
1 |3 6
2 |4
3 |5
4 |6
5 |7 2

julia> selct(t, a=d->mean(d[:a]), b=d->sum(d[:b] .* 2), by=Any[:d=>d->d[:b] .* 2], by=[:c])
4 x 3 LabeledArray

c |x     |y     |z
--+------+------+-------
d |a   b |a   b |a    b
--+------+------+-------
  |4.0 0 |5.5 0 |
2 |1.0 2 |7.0 2 |8.0  2
4 |2.0 4 |      |9.0  4
6 |3.0 6 |      |10.0 6
```

"""
function selct end

selct(t::LabeledArray, args...; kwargs...) = selct_untyped(t, args...; kwargs...)
selct(t::DictArray, args...; kwargs...) = selct_untyped(t, args...; kwargs...)
selct_untyped(t, args...; kwargs...) = begin
  a = Any[args...]
  b = []
  c = []
  for (k,v) in kwargs
    if k == :by
      push!(b, isa(v, AbstractArray) ? v : [v])
    elseif k == :where
      push!(c, (isa(v, AbstractArray) ? v : [v])...)
    else
      push!(a, k=>v)
    end
  end
  sel(selectfunc, t; c=c, b=b, a=a)
end

"""

`@update(t, args...)`

Similar to `select` macro, but is used to update and create a new `LabeledArray` or `DictArray` from the original one.
The main difference from the `select` macro is that it keeps the original fields intact, unless
directed otherwise, whereas the `select` macro only chooses the fields that are explicitly specified.

##### Arguments
Below `t'` is an object such that `t'[k]` for a field name k gives the corresponding array
for the field `k` in the table `t` at the coordinates selected so far.

* `t` : a `LabeledArray` or `DictArray`.
* `args...` : each argument can be field names, fieldname=>(`t'`->nullable array function) pair or `where[...]`. A fieldname=function can be used instead of the pair notation if the fieldname is a symbol. `by[...]` has (`t'`->nullable array) as its elements and used as in grouping when updateing. Multiple `by[...]` are simply combined. `where[...]` has
(`t'`->nullable boolean array) functions inside and chooses the appropriate portion of the original array. Multiple `where[...]` will simply be combined.

##### Function Specification

Note that a function (`t'`->nullable array) is expressed by some **expression** with variable names with underscores. The **expression** is converted into `t''->**expression**. Symbols with underscores are converted in the following way:

* `_k` : translates to `t'[k]`. The field name `k` should be a symbol.
* `_!k` : translates to `isna(t')[k]`. It gives a boolean array to denote whether an element is null. The field name `k` should be a symbol.
* `_` : translates to `t'` if none of the above is applicable. It is useful when the field name is not a symbol.

##### Return

An updated array of the same type as `t`.

##### Examples

```julia
julia> t = @larr(a=1:10, b=[1,2,3,NA,NA,NA,1,1,2,3], c=[:x,:x,:x,:x,:y,:y,:y,:z,:z,:z])
10 LabeledArray

   |a  b c
---+-------
1  |1  1 x
2  |2  2 x
3  |3  3 x
4  |4    x
5  |5    y
6  |6    y
7  |7  1 y
8  |8  1 z
9  |9  2 z
10 |10 3 z

julia> @update(t, a=_a .+ 100, d=_a .* _b)
10 LabeledArray

   |a   b c d
---+-----------
1  |101 1 x 1
2  |102 2 x 4
3  |103 3 x 9
4  |104   x
5  |105   y
6  |106   y
7  |107 1 y 7
8  |108 1 z 8
9  |109 2 z 18
10 |110 3 z 30

julia> @update(t, a=_a .+ 100, d=_a .* _b, where[~isna(_b)])
10 LabeledArray

   |a   b c d
---+-----------
1  |101 1 x 1
2  |102 2 x 4
3  |103 3 x 9
4  |4     x
5  |5     y
6  |6     y
7  |107 1 y 7
8  |108 1 z 8
9  |109 2 z 18
10 |110 3 z 30

julia> @update(t, a=sum(_a), d=reverse(_a .* _b), where[~isna(_b)], by[:b])
10 LabeledArray

   |a  b c d
---+----------
1  |16 1 x 8
2  |11 2 x 18
3  |13 3 x 30
4  |4    x
5  |5    y
6  |6    y
7  |16 1 y 7
8  |16 1 z 1
9  |11 2 z 4
10 |13 3 z 9
```

"""
macro update(t, args...)
  select_inner(:updatefunc, t, args...)
end

"""

`update(t, agg... [; by=[...]..., where=[...]...])`

Similar to `selct`, but is used to update and create a new `LabeledArray` or `DictArray` from the original one.
The main difference from the `selct` function is that it keeps the original fields intact, unless
directed otherwise, whereas the `select` macro only chooses the fields that are explicitly specified.

##### Arguments

Below `t'` is an object such that `t'[k]` for a field name `k` gives the corresponding array
for the field `k` in the table `t` at the coordinates selected so far.

* `t` : a LabeledArray or `DictArray`.
* `agg...` : each argument can be field names, fieldname=>(t'->nullable array function) pair, `by[...]` or `where[...]`. A fieldname=function can be used instead of the pair notation if the fieldname is a symbol.
* `by=[...]`: `by[...]` has (`t'`->nullable array) as its elements and used as in grouping when updateing. Multiple `by[...]` are simply combined.
* `where=[...]`: has (`t'`->nullable boolean array) functions inside and chooses the appropriate portion of the original array. Multiple `where=[...]` will simply be combined.

##### Return

An updated array of the same type as `t`.

##### Examples

```julia
julia> t = @larr(a=1:10, b=[1,2,3,NA,NA,NA,1,1,2,3], c=[:x,:x,:x,:x,:y,:y,:y,:z,:z,:z])
10 LabeledArray

   |a  b c
---+-------
1  |1  1 x
2  |2  2 x
3  |3  3 x
4  |4    x
5  |5    y
6  |6    y
7  |7  1 y
8  |8  1 z
9  |9  2 z
10 |10 3 z

julia> update(t, a=d->d[:a] .+ 100, d=d->d[:a] .* d[:b])
10 LabeledArray

   |a   b c d
---+-----------
1  |101 1 x 1
2  |102 2 x 4
3  |103 3 x 9
4  |104   x
5  |105   y
6  |106   y
7  |107 1 y 7
8  |108 1 z 8
9  |109 2 z 18
10 |110 3 z 30

julia> update(t, a=d->d[:a] .+ 100, d=d->d[:a] .* d[:b], where=[d-> ~isna(d[:b])])
10 LabeledArray

   |a   b c d
---+-----------
1  |101 1 x 1
2  |102 2 x 4
3  |103 3 x 9
4  |4     x
5  |5     y
6  |6     y
7  |107 1 y 7
8  |108 1 z 8
9  |109 2 z 18
10 |110 3 z 30

julia> update(t, a=d->sum(d[:a]), d=d->reverse(d[:a] .* d[:b]), where=[d-> ~isna(d[:b])], by=[:b])
10 LabeledArray

   |a  b c d
---+----------
1  |16 1 x 8
2  |11 2 x 18
3  |13 3 x 30
4  |4    x
5  |5    y
6  |6    y
7  |16 1 y 7
8  |16 1 z 1
9  |11 2 z 4
10 |13 3 z 9
```

"""
function update end

update(t::LabeledArray, args...; kwargs...) = update_untyped(t, args...; kwargs...)
update(t::DictArray, args...; kwargs...) = update_untyped(t, args...; kwargs...)
update_untyped(t, args...; kwargs...) = begin
  a = Any[args...]
  b = []
  c = []
  for (k,v) in kwargs
    if k == :by
      push!(b, isa(v, AbstractArray) ? v : [v])
    elseif k == :where
      push!(c, (isa(v, AbstractArray) ? v : [v])...)
    else
      push!(a, k=>v)
    end
  end
  sel(updatefunc, t; c=c, b=b, a=a)
end

# a helper function to decipher expressions for the select and update macros.
select_inner(func, t, args...) = begin
  cexprs = Any[]
  bexprs = Any[]
  aexprs = Any[]
  for arg in args
    if :head in fieldnames(arg)
      if arg.head == :kw || arg.head == :(=>) || arg.head == :(=)
        push!(aexprs, replace_expr(arg))
      elseif arg.head == :ref && arg.args[1] == :where
        push!(cexprs, map(replace_expr, arg.args[2:end])...)
      elseif arg.head == :ref && arg.args[1] == :by
        push!(bexprs, map(x -> replace_expr(lift_to_dict(x)), arg.args[2:end]))
      else
        # when arg is a quoted symbol.
        push!(aexprs, replace_expr(lift_to_dict(arg)))
      end
    else
      push!(aexprs, replace_expr(lift_to_dict(arg)))
    end
  end
  cexprvec = if isempty(cexprs); :([]) else Expr(:vect, cexprs...) end
  bexprvec = if isempty(bexprs); :([]) else Expr(:ref, :Any, map(v->Expr(:vect, v...), bexprs)...) end
  aexprvec = if isempty(aexprs); :([]) else Expr(:vect, aexprs...) end
  quote
    sel($func, $(esc(t)), c=$cexprvec, b=$bexprvec, a=$aexprvec)
  end
end

const verbatim_magic_symbol = symbol("__verbatim__")

lift_to_dict(expr) = begin
  if :head in fieldnames(expr) &&
     (expr.head == :kw || expr.head == :(=>) || expr.head == :(=))
    expr
  else
    Expr(:(=>), expr, Expr(:ref, verbatim_magic_symbol, expr))
  end
end

"""

`replace_expr(expr)`

Create a function expression from a domain expression.

##### Expressions
Below t' is an object such that t'[k] for a field name k gives the corresponding array
for the field k in the table t at the coordinates selected so far.

* `_k` : translates to `t'[k]`. The field name `k` should be a symbol.
* `__k` : translates to `igna(t')[k]`. It ignores the null elements. The null elements are replaced with arbitrary values, so make sure there is no null value in the array if you want to use it. The field name `k` should be a symbol.
* `_!k` : translates to `isna(t')[k]`. It gives a boolean array to denote whether an element is null. The field name `k` should be a symbol.
* `_` : translates to `t'` if none of the above is applicable. It is useful when the field name is not a symbol.

"""
replace_expr(expr) = begin
  noescape_symbols = [:none, :LineNumberNode, :QuoteNode, :line, :quote]
  result = if :head in fieldnames(expr) && expr.head == :kw || expr.head == :(=>) || expr.head == :(=)
    # this is for a or b case.
    key = if expr.head == :(=) || expr.head == :kw; quote_symbol(expr.args[1]) else expr.args[1] end #quote_symbol(expr.args[1])
    if :head in fieldnames(expr.args[2]) && expr.args[2].head == :ref && expr.args[2].args[1] == verbatim_magic_symbol
      Expr(:(=>), esc(key), esc(expr.args[2].args[2]))
    else
      value = replace_expr_inner(expr.args[2], noescape_symbols...)
      Expr(:(=>), esc(key), Expr(:(->), :d, Expr(:call, :nalift, value)))
    end
  else
    # this is for c case.
    Expr(:(->), :d, replace_expr_inner(expr, noescape_symbols...))
  end
  result
end

replace_expr_inner(expr, noescape_symbols...) = begin
  strexpr = string(expr)
  if isa(expr, Symbol) && startswith(strexpr, "__")
    Expr(:ref, Expr(:call, :igna, :d), QuoteNode(symbol(strexpr[3:endof(strexpr)])))
  elseif isa(expr, Symbol) && startswith(strexpr, "_!")
    Expr(:ref, Expr(:call, :isna, :d), QuoteNode(symbol(strexpr[3:endof(strexpr)])))
  elseif isa(expr, Symbol) && strexpr == "_"
    :d
  elseif isa(expr, Symbol) && startswith(strexpr, '_')
    Expr(:ref, :d, QuoteNode(symbol(strexpr[2:endof(strexpr)])))
  elseif :head in fieldnames(expr) && expr.head==:macrocall && expr.args[1] == symbol("@rap")
    # if it is @rap, give precedence to it.
    replace_expr_inner(recursive_rap(expr.args[2:end]...), noescape_symbols...)
  elseif :head in fieldnames(expr) && expr.head==:ref && expr.args[1]==symbol("__")
    Expr(:ref, Expr(:call, :igna, :d), esc(expr.args[2]))
  elseif :head in fieldnames(expr) && expr.head==:ref && expr.args[1]==symbol("_!")
    Expr(:ref, Expr(:call, :isna, :d), esc(expr.args[2]))
  elseif :head in fieldnames(expr) && expr.head==:ref && expr.args[1]==:_
    #TODO need to fix this so that it works with an array of field names.
    Expr(:ref, replace_expr_inner(expr.args[1], noescape_symbols...), map(x->isa(x,QuoteNode) || (:head in fieldnames(x) && x.head==:quote) ? x : replace_expr_inner(x, noescape_symbols...), expr.args[2:end])...)
    #Expr(:ref, expr.args[1], map(x->isa(x, Symbol) ? QuoteNode(x) : replace_expr_inner(x), expr.args[2:end])...)
  elseif :head in fieldnames(expr) && expr.head == :(->)
    noescsymbs = get_function_arguments_exprs(expr.args[1])
    Expr(:(->), expr.args[1], replace_expr_inner(expr.args[2], noescsymbs..., noescape_symbols...))
  else
    if isa(expr, QuoteNode)
      expr
    elseif :head in fieldnames(expr) && expr.head in noescape_symbols
      expr
    elseif :args in fieldnames(expr)
      newargs = map(x->replace_expr_inner(x, noescape_symbols...), expr.args)
      expr.args = newargs
      expr
    elseif expr in noescape_symbols
      expr
    else
      esc(expr)
    end
  end
end

get_function_arguments_exprs(expr) = begin
  if :head in fieldnames(expr) && expr.head == :tuple
    expr.args
  else
    (expr,)
  end
end

# N is the original array dimensions.
update_by_array{N}(t::AbstractArray{TypeVar(:T),N}, a, indices::Vector{NTuple{N,Int}}, byvec, bymap) = begin
  label = bymap[1]
  # assumes there is at least one label for each axis.
  reslen = length(bymap[1][1])
  indexvec::Vector{Int} = bymap[2]
  indcountvec::Vector{Int} = fill(0, reslen)
  for i in eachindex(indices)
    ids = indexvec[i]
    indcountvec[ids] += 1
  end
  subarrays_buf = Array(NTuple{N,Int}, length(indices))
  index_subarrays = similar(indcountvec, SubArray{NTuple{N,Int},1,Array{NTuple{N,Int},1},Tuple{UnitRange{Int}},1})
  indsumsofar::Int = 0
  lenindcountvec::Int = length(indcountvec)
  for i in 1:lenindcountvec
    prevsum = indsumsofar
    @inbounds indsumsofar += indcountvec[i]
    @inbounds index_subarrays[i] = sub(subarrays_buf, prevsum+1:indsumsofar)
  end
  # change the meaning of indcountmat, not to create another buffer.
  fill!(indcountvec, 0)
  counter::Int = 1
  for i in eachindex(indices)
    ids = indexvec[i]
    current_count = indcountvec[ids] + 1
    indcountvec[ids] = current_count
    subarray = index_subarrays[ids]
    subarray[current_count] = indices[i]
  end
  result = copy(t.data.data)
  for (k,onea) in a
    trial_indices = index_subarrays[1]
    trial_aggvec = onea(t, trial_indices)
    result = create_field_for_update!(result, k, trial_aggvec)
    update_agg!(result, k, trial_aggvec, trial_indices)
    start_index = 2

    for i in start_index:length(index_subarrays)
      indices = index_subarrays[i]
    #for indices in index_subarrays
      aggvec = onea(t, indices)
      update_agg!(result, k, aggvec, indices)
    end
  end
  LabeledArray(DictArray(result), t.axes)
end

create_field_for_update!{T}(result::LDict, k, ::AbstractArray{Nullable{T}}) = begin
  k_in_result = haskey(result, k)
  repval = values(result)[1]
  if !k_in_result
    v = similar(repval, Nullable{T}, size(repval))
    setna!(v)
    merge(result, LDict(k=>v))
  elseif !isa(T, eltype(eltype(result[k])))
    v = similar(repval, Nullable{promote_type(T,eltype(eltype(result[k])))}, size(repval))
    copy!(v, result[k])
    merge(result, LDict(k=>v))
  else
    result
  end
end
create_field_for_update!{T}(result::LDict, k, ::Nullable{T}) = begin
  k_in_result = haskey(result, k)
  repval = values(result)[1]
  if !k_in_result
    v = similar(repval, Nullable{T}, size(repval))
    setna!(v)
    merge(result, LDict(k=>v))
  elseif !isa(T, eltype(eltype(result[k])))
    v = similar(repval, Nullable{promote_type(T,eltype(eltype(result[k])))}, size(repval))
    copy!(v, result[k])
    merge(result, LDict(k=>v))
  else
    result
  end
end
create_field_for_update!{T}(result::LDict, k, ::T) = create_field_for_update!(result, k, Nullable{Any})

update_agg!{N}(result::LDict, k, aggvec::AbstractVector, indices::AbstractVector{NTuple{N,Int}}) = begin
  @assert(length(aggvec) == length(indices))
  i = 1
  resultk = result[k]
  for coords in indices
    resultk[coords...] = aggvec[i]
    i += 1
  end
end

update_agg!{N}(result::LDict, k, aggnum, indices::AbstractVector{NTuple{N,Int}}) = update_agg!(result, k, Nullable(aggnum), indices)
update_agg!{N,V}(result::LDict, k, aggnum::Nullable{V}, indices::AbstractVector{NTuple{N,Int}}) = begin
  i = 1
  resultk = result[k]
  for coords in indices
    resultk[coords...] = aggnum
    i += 1
  end
end


text:
# The purpose of this script is to compare the filtering capability of the single snapshots
# estimates and the gas filter on a given

using StatsBase, JLD,  Hwloc,  DynNets, StaticNets, StatsFuns,Utilities

# static Estimates for increasing N ---------------------------------------------

#need to choose how to distribute the unconditional means. I want to stay
# in the case of chaterjee Diaconis Sly, hence have fixed maximum value for the
# parameters that does not increase with N
fractionParTV = 0.05;

Nrange = 50:100:800
Nsample = 50
degMin = 5
degMax = 30

midInd(vec::Array{<:Real,1}) = round(Int,length(vec)/2)

save_fold = "./data/estimatesTest/asympTest/"

data = load(save_fold*"asympTestStaticLarge.jld")
estFits = data["estFits"]
dgpFits = data["dgpFits"]
N,Scal = size(dgpFits)
Sam = size(estFits[1,1])[2]
expDegsDgp = copy(dgpFits)
for n=1:N,scal=1:Scal expDegsDgp[n,scal] = expValStats(fooErgmDirBin1,dgpFits[n,scal]) end
expDegsEst = copy(estFits)
for n=1:N,scal=1:Scal,sam =1:Sam  expDegsEst[n,scal][:,sam] = expValStats(fooErgmDirBin1,estFits[n,scal][:,sam]) end

# plot the Degrees in the dgp
for n=1:N
    scal = 3
    plot(expDegsDgp[n,scal][1:midInd(expDegsDgp[n,scal])])
 end
 grid()

# plot densities
for scal=[1,3]
     dens = [sum(expDegsDgp[n,scal])/(length(expDegsDgp[n,scal])*(length(expDegsDgp[n,scal])-1)) for n=1:N]
     size = [length(expDegsDgp[n,scal]) for n=1:N]
     plot(size,dens ,".-")
  end
  grid()
  xlabel("number of nodes")
  ylabel("density")

for n=1:N
    scal = 3

 diffDegs = (expDegsDgp[n,scal] .- expDegsEst[n,scal])./expDegsDgp[n,scal]
  sampMenDiff = meanSq(diffDegs,2)
 plot(expDegsDgp[n,scal][1:midInd(expDegsDgp[n,scal])],n*0.2 + sampMenDiff[1:midInd(sampMenDiff)],".")
 end
 yticks(0.2:0.2:N*0.2)
 grid()

for n=1:N
     scal = 3
  diffDegs = (expDegsDgp[n,scal] .- expDegsEst[n,scal])./expDegsDgp[n,scal]
  sampStdDiff = squeeze(std(diffDegs,2),2)
  plot(expDegsDgp[n,scal][1:midInd(sampStdDiff)],n*0.2 + sampStdDiff[1:midInd(sampStdDiff)],".")
 end
  yticks(0.2:0.2:N*0.2)
  grid()

##

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
# ev = CPEG.CPEGWorkspace()
#
# Rm = ev.params.gravity.R
# r0 = SA[Rm+125e3, 0.0, 0.0] #Atmospheric interface at 125 km altitude
# V0 = 5.845e3 #Mars-relative velocity at interface of 5.845 km/sec
# γ0 = -15.474*(pi/180.0) #Flight path angle at interface
# v0 = V0*SA[sin(γ0), cos(γ0), 0.0]
# σ0 = deg2rad(90)
#
# r0sc,v0sc = CPEG.scale_rv(ev.scale,r0,v0)
#
# dt = 1.0/ev.scale.tscale
#
# x0 = SA[r0sc[1],r0sc[2],r0sc[3],v0sc[1],v0sc[2],v0sc[3],σ0]
#
# N = 100
# U = [SA[0.0] for i = 1:N-1]
# #
# X1,U1 = CPEG.rollout(ev, x0, U, dt)
#
# @test length(X1) == 113
# @test (norm((X1[end])[1:3]) * ev.scale.dscale - ev.params.gravity.R) / 1000.0 < 10.0
# @test (norm((X1[end-1])[1:3]) * ev.scale.dscale - ev.params.gravity.R) / 1000.0 > 10.0
#
#
# # now update the scaling
#
# ev.scale.dscale = 1.0
# ev.scale.tscale = 1.0
# ev.scale.uscale = 1.0
#
# x0 = SA[r0[1],r0[2],r0[3],v0[1],v0[2],v0[3],σ0]
# dt = 1.0/ev.scale.tscale
# X,U = CPEG.rollout(ev, x0, U, dt)
#
# @test length(X) == 113
# @test (norm((X[end])[1:3]) * ev.scale.dscale - ev.params.gravity.R) / 1000.0 < 10.0
# @test (norm((X[end-1])[1:3]) * ev.scale.dscale - ev.params.gravity.R) / 1000.0 > 10.0

# # test scaling
# for i = 1:length(X)
#
#     r1 = SVector{3}(X1[i][1:3])

function test_rollout(ev,r0,v0,σ0)

    r0sc,v0sc = CPEG.scale_rv(ev.scale,r0,v0)

    # dt = 1.0/ev.scale.tscale
    ev.dt = 1.0

    x0 = SA[r0sc[1],r0sc[2],r0sc[3],v0sc[1],v0sc[2],v0sc[3],σ0]

    X, U = CPEG.rollout(ev, x0, ev.U)

    # @show length(X)

    @test length(X) == 241
    @test length(U) == 241 - 1
    @test (norm((X[end])[1:3]) * ev.scale.dscale - ev.params.gravity.R)  < 10.0e3
    @test (norm((X[end-1])[1:3]) * ev.scale.dscale - ev.params.gravity.R) > 10.0e3

    return X
end

ev1 = CPEG.CPEGWorkspace()

Rm = ev1.params.gravity.R
r0 = SA[Rm+125e3, 0.0, 0.0] #Atmospheric interface at 125 km altitude
V0 = 5.845e3 #Mars-relative velocity at interface of 5.845 km/sec
γ0 = -15.474*(pi/180.0) #Flight path angle at interface
v0 = V0*SA[sin(γ0), cos(γ0), 0.0]
σ0 = deg2rad(42)

X1 = test_rollout(ev1,r0,v0,σ0)

ev2 = CPEG.CPEGWorkspace()

ev2.scale.dscale = 1.0
ev2.scale.tscale = 1.0
ev2.scale.uscale = 1.0

X2 = test_rollout(ev2,r0,v0,σ0)

@test length(X1) == length(X2)
for i = 1:length(X1)

    r1_sc = SVector{3}(X1[i][1:3])
    v1_sc = SVector{3}(X1[i][4:6])

    r1,v1 = CPEG.unscale_rv(ev1.scale,r1_sc,v1_sc)

    r2_sc = SVector{3}(X2[i][1:3])
    v2_sc = SVector{3}(X2[i][4:6])

    r2,v2 = CPEG.unscale_rv(ev2.scale,r2_sc,v2_sc)

    @test r1 ≈ r2 rtol = 1e-6
    @test v1 ≈ v2 rtol = 1e-6
end