text:
"""
Parse `var::T` or `var`
"""
function varT(input)

    # Type Checker
    if input isa Symbol #  var

        var = input

        args_type = :Any
        args_field =  var
        args_param = var

        exprargsI = var # remove = defVal for type definition

    elseif input isa Expr && input.head == :(::) && input.args[1] isa Symbol # var::T

        var = input.args[1]
        varType = input.args[2] # Type

        args_type = varType
        args_field =  var
        args_param = var

        exprargsI = input # remove =defVal for type definition

    else
        # something else, e.g. inline inner constructor
        #   F(...) = ...
        args_type, args_field, args_param, exprargsI = missing, missing, missing, nothing
    end

    return args_type, args_field, args_param, exprargsI
end

"""
Parse `var::T = defVal` or `var = defVal`
"""
function varTdefVal(input)
    input.head == :(kw) || input.head == :(=) || throw(ArgumentError("Input should be in form of `var/var::T = defVal`. You passed $input"))

    lhs = input.args[1] # var or var::T
    args_type, args_field, args_param, exprargsI = varT(lhs)

    args_defaultvalue = input.args[2] # defVal
    args_param = Expr(:kw, args_field, args_defaultvalue)

    return args_defaultvalue, args_type, args_param, args_field, exprargsI
end

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

"""
VarStruct definer function. It also serves as a parser.
# Examples
```julia
varstruct_call_parse(
    Main,
    quote
        Person(
            name::String,
            number::Float64 = 20.0,
        )
    end
)
```
"""
function varstruct_call_parse(modul::Module, expr::Expr)
    # The following code is taken and simplified from the `aml_parse` function of AcuteML.jl (@aminya):

    if expr.head == :block
        expr = expr.args[2]
    end
    if expr.head !== :call
        error("Invalid usage of `@var`")
    end

    T = expr.args[1] # Type name

    if !isa(T, Symbol)
        error("Invalid name for the type. Change `$T` to a Symbol")
    end

    data = expr.args[2:end] # arguments of the type
    argsnum = length(data)

    args_param = Vector{Union{Expr,Symbol}}(undef, argsnum) # Expr(:parameters)[]
    args_field = Vector{Union{Symbol, String}}(undef, argsnum)
    args_defaultvalue = Vector{Any}(undef, argsnum)
    fill!(args_defaultvalue, Unset())
    args_type = Vector{Union{Type, Symbol, Expr}}(undef, argsnum)

    for iArg = 1:argsnum # iterating over arguments of each type argument

        ei = data[iArg] # type argument element i
        iExprArgs = iArg + 1
        ################################################################
        # Arguments
        ################################################################
        # No defVal
        if (ei isa Symbol) || (ei.head == :(::)) # var/var::T

            error("You should provide an initial value for $ei")

            args_type[iArg], args_field[iArg], args_param[iArg], exprargsI = varT(ei)
            !isnothing(exprargsI) ? expr.args[iExprArgs] = exprargsI : nothing
        ################################################################
        # Function not provided
        elseif ei.head == :(kw) # def value provided
             # var/var::T = defVal

             args_defaultvalue[iArg], args_type[iArg], args_param[iArg], args_field[iArg], exprargsI = varTdefVal(ei)
             !isnothing(exprargsI) ? expr.args[iExprArgs] = exprargsI : nothing
        ################################################################
        # Needs reevaluation
        elseif ei.head == :call
            error("Chaining not supported yet")
            expr.args[iExprArgs],  args_field[iArg], args_defaultvalue[iArg], args_type[iArg],  args_param[iArg], T  = varstruct_call_parse(modul, expr)
            # should call creator here:

        elseif ei.head == :block
            error("Chaining not supported yet")
            expr.args[iExprArgs],  args_field[iArg], args_defaultvalue[iArg], args_type[iArg],  args_param[iArg], T  = varstruct_call_parse(modul, expr.args[2])
            # should call creator here:
        else
            error("Syntax is not supported")
        end  # end ifs
    end # endfor

    return expr, args_field, args_defaultvalue, args_type, args_param, T
end

################################################################
"""
VarStruct declaration function. Use this syntax when you only need to declare the struct and not fully initialize it.
# Examples
```julia
varstruct_struct_parse(
    Main,
    quote
        struct Person
            name::String
            number::Float64 = 20.0
        end
    end
)
```
"""
function varstruct_struct_parse(modul::Module, expr::Expr)
    # The following code is taken and simplified from the aml_parse function of AcuteML.jl (@aminya):

    if expr.head == :block
        expr = expr.args[2]
    end
    if expr.head !== :struct
        error("Invalid usage of `@var`")
    end

    is_struct_mutable = expr.args[1] # for future use
    T = expr.args[2] # Type name +(curly braces)

    if !isa(T, Symbol)
        error("Invalid name for the type. Change `$T` to a Symbol")
    end

    argsexpr = expr.args[3] # arguments of the type

    areargs_inds = findall(x->!(isa(x, LineNumberNode)), argsexpr.args)
    data = argsexpr.args[areargs_inds]
    argsnum = length(areargs_inds)

    args_param = Vector{Union{Expr,Symbol}}(undef, argsnum) # Expr(:parameters)[]
    args_field = Vector{Union{Symbol, String}}(undef, argsnum)
    args_defaultvalue = Vector{Any}(undef, argsnum)
    fill!(args_defaultvalue, Unset())
    args_type = Vector{Union{Type, Symbol, Expr}}(undef, argsnum)
    args_check = Vector{Union{Missing, Symbol, Function}}(missing, argsnum)
    # struct_check = Array{Union{Missing, Symbol, Function},0}(missing)

    for iArg = 1:argsnum # iterating over arguments of each type argument

        i = areargs_inds[iArg] # used for argsexpr.args[i]

        ei = data[iArg] # type argument element i

        ################################################################
        # Arguments
        # No defVal
        if (ei isa Symbol) || (ei.head == :(::)) # var/var::T

            args_type[iArg], args_field[iArg], args_param[iArg], exprargsI = varT(ei)
            !isnothing(exprargsI) ? argsexpr.args[i] = exprargsI : nothing

        ########################
        # Function provided
        elseif  ei.head == :tuple

            # var/var::T = defVal, F
            if isa(ei.args[2], Union{Function, Symbol})  # function provided

                # var/var::T = defVal
                args_defaultvalue[iArg], args_type[iArg], args_param[iArg], args_field[iArg], exprargsI = varTdefVal(ei.args[1])
                !ismissing(exprargsI) ? argsexpr.args[i] = exprargsI : nothing

                 # F
                args_check[iArg] = ei.args[2]
            else
                error("Unkown syntax.")
            end

        ################################################################
        # Function not provided
        elseif ei.head == :(=) # def value provided
             # var/var::T = defVal

             args_defaultvalue[iArg], args_type[iArg], args_param[iArg], args_field[iArg], exprargsI = varTdefVal(ei)
             !isnothing(exprargsI) ? argsexpr.args[i] = exprargsI : nothing

        # Needs reevaluation
        elseif ei.head == :struct
             error("Chaining not supported yet")
             # TODO: only works for structs
             argsexpr.args[i], args_field[iArg], args_defaultvalue[iArg], args_type[iArg], args_param[iArg], args_check[iArg], _, _  = varstruct_struct_parse(modul, expr)

        ################################################################
        else
            error("Syntax is not supported")
        end # end ifs
    end # endfor

    return argsexpr, args_field, args_defaultvalue, args_type, args_param, args_check, is_struct_mutable, T
end


text:
# Join node.

mutable struct JoinNode <: TabularNode
    over::Union{SQLNode, Nothing}
    joinee::SQLNode
    on::SQLNode
    left::Bool
    right::Bool

    JoinNode(; over = nothing, joinee, on, left = false, right = false) =
        new(over, joinee, on, left, right)
end

JoinNode(joinee; over = nothing, on, left = false, right = false) =
    JoinNode(over = over, joinee = joinee, on = on, left = left, right = right)

JoinNode(joinee, on; over = nothing, left = false, right = false) =
    JoinNode(over = over, joinee = joinee, on = on, left = left, right = right)

"""
    Join(; over = nothing, joinee, on, left = false, right = false)
    Join(joinee; over = nothing, on, left = false, right = false)
    Join(joinee, on; over = nothing, left = false, right = false)

`Join` correlates two input datasets.

```sql
SELECT ...
FROM \$over
JOIN \$joinee ON \$on
```

# Examples

```jldoctest
julia> person = SQLTable(:person, columns = [:person_id, :location_id]);

julia> location = SQLTable(:location, columns = [:location_id, :state]);

julia> q = person |>
           Join(:location => location,
                Get.location_id .== Get.location.location_id) |>
           Select(Get.person_id, Get.location.state);

julia> print(render(q))
SELECT "person_1"."person_id", "location_1"."state"
FROM "person" AS "person_1"
JOIN "location" AS "location_1" ON ("person_1"."location_id" = "location_1"."location_id")
```
"""
Join(args...; kws...) =
    JoinNode(args...; kws...) |> SQLNode

"""
An alias for `Join(...; ..., left = true)`.
"""
LeftJoin(args...; kws...) =
    Join(args...; kws..., left = true)

dissect(scr::Symbol, ::typeof(Join), pats::Vector{Any}) =
    dissect(scr, JoinNode, pats)

function PrettyPrinting.quoteof(n::JoinNode, ctx::QuoteContext)
    ex = Expr(:call, nameof(Join))
    if !ctx.limit
        push!(ex.args, quoteof(n.joinee, ctx))
        push!(ex.args, quoteof(n.on, ctx))
        if n.left
            push!(ex.args, Expr(:kw, :left, n.left))
        end
        if n.right
            push!(ex.args, Expr(:kw, :right, n.right))
        end
    else
        push!(ex.args, :…)
    end
    if n.over !== nothing
        ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
    end
    ex
end

label(n::JoinNode) =
    label(n.over)

rebase(n::JoinNode, n′) =
    JoinNode(over = rebase(n.over, n′),
             joinee = n.joinee, on = n.on, left = n.left, right = n.right)

ext: jl


text:
abstract type StaticInjection <: Device end
supports_services(d::Device) = true


text:
"""
Initialize an model with zeros, given mgrid and names.
"""
function Base.zero(::Type{Medium}, mgrid, names = [:vp, :rho])
    nvec = length.(mgrid)
    return Medium(
        map( StepRangeLen{Float64}, mgrid),
        # broadcast(
            # x -> convert(
            #    ,
                # x,
            # ),
            # (mgrid),
        # ),
        NamedArray([zeros(Float64, nvec...) for i in names], (names,)),
        NamedArray([fill(0.0, 2) for i in names], (names,)),
        NamedArray([0.0 for i in names], (names,)),
    )
end
Medium(mgrid, names = [:vp, :rho]) = zero(Medium, mgrid, names)
Medium(mgrid, m::Medium) = zero(Medium, mgrid, names(m.m)[1])

"""
Return true if a `Medium` model is just allocated with zeros.
"""
function Base.iszero(mod::Medium)
    if (any(broadcast(iszero, mod.bounds)) && any(broadcast(iszero, mod.ref)))
        return true
    else
        # vp and rho cannot be zero
        for field in [:vp, :rho]
            if (field ∈ names(mod.m)[1])
                if (any(iszero.(mod[field])))
                    return true
                end
            end
        end
    end
    return false
end

"""
Return true if a `Medium` model is just allocated with zeros.
"""
function Base.isassigned(mod::Medium, fields = [:vp, :rho, :vs])
    for field in fields
        @assert field ∈ [:vp, :rho, :vs]
    end

    return any(broadcast(iszero, mod.bounds)) && any(broadcast(iszero, mod.ref))
end

"""
Return a similar model (size, bounds and reference values) to the input model, used for allocation.
"""
function Base.similar(mod::Medium)
    mod_new = deepcopy(mod)
    fill!(mod_new)
    return mod_new
end

"Compare if two `Medium` models are equal"
function Base.isequal(mod1::T, mod2::T) where {T<:Union{Medium}}
    fnames = fieldnames(typeof(mod1))
    vec = [(isequal(getfield(mod1, name), getfield(mod2, name))) for name in fnames]
    return all(vec)
end

"""
Return if two `Medium` models have same dimensions and bounds.
"""
function Base.size(mod::Medium)
    return length.(mod.mgrid)
end

"""
Return the number of dimensions of `Medium`
"""
function Base.ndims(mod::Medium)
    return length(mod.mgrid)
end

"""
Copy for `Medium`. The media should have same bounds and sizes.
Doesn't allocate any memory. Note that only [:vp, :rho, :vs] fields are stored in structs.
"""
function Base.copyto!(modo::Medium, mod::Medium)
    for name in [:vp, :vs, :rho]
        ((name ∈ names(modo.m)[1]) & (name ∈ names(mod.m)[1])) &&
            (copyto!(modo.m[name], mod.m[name]))
    end
    return modo
end

"""
Fill medium with reference medium parameters to make it homogeneous.
```julia
fill!(medium)
```
"""
function Base.fill!(mod::Medium)
    for s in names(mod.m)[1]
        fill!(mod.m[s], mod.ref[s])
    end
    return mod
end

"""
Print information about `Medium`
"""
function Base.show(io::Base.IO, mod::Medium)
    println(
        io,
        "> ",
        ndims(mod),
        "-D Medium with ",
        dim_names(ndims(mod)),
        " in ",
        [[m[1], m[end]] for m in mod.mgrid],
    )
    println(io, "> number of samples:\t", length.(mod.mgrid))
    println(io, "> sampling intervals:\t", step.(mod.mgrid))
    try
        println(io, "> vp:\t", "min\t", minimum(mod[:vp]), "\tmax\t", maximum(mod[:vp]))
    catch
    end
    try
        println(io, "> vs:\t", "min\t", minimum(mod[:vs]), "\tmax\t", maximum(mod[:vs]))
    catch
    end
    try
        println(io, "> rho:\t", "min\t", minimum(mod[:rho]), "\tmax\t", maximum(mod[:rho]))
    catch
    end
    println(io, "> Bounds:")
    println(io, mod.bounds)
end

"""
Re-parameterization routine
that modifies the fields
`χvp` and `χrho` of an input seismic model
using two input vectors.

# Arguments

* `mod::Medium` : to be updated
* `x1::Array{Float64,2}` : contrast of inverse bulk modulus
* `x2::Array{Float64,2}` : contrast of inverse density
* `fields:::Vector{Symbol}` : [:χKI, :χrhoI]
"""
function Base.copyto!(
    mod::Medium,
    x::AbstractVector,
    fields::Vector{Symbol} = [:χKI, :χrhoI, :null],
)
    iszero(mod) ? error("mod cannot be zero") : nothing
    nznx = prod(length.(mod.mgrid))
    (length(x) ≠ (count(fields .≠ :null) * nznx)) && error("size x")
    if (fields == [:χKI, :χrhoI, :null])
        @inbounds for i = 1:nznx
            K = inv(χ(x[i], mod.ref[:KI], -1))
            rho = inv(χ(x[nznx+i], mod.ref[:rhoI], -1))
            mod.m[:vp][i] = sqrt(K * inv(rho))
            mod.m[:rho][i] = rho
        end
    elseif (fields == [:χKI, :null, :null])
        @inbounds for i = 1:nznx
            K = inv(χ(x[i], mod.ref[:KI], -1))
            rho = mod.m[:rho][i]
            mod.m[:vp][i] = sqrt(K * inv(rho))
            mod.m[:rho][i] = rho
        end
    elseif (fields == [:χvp, :χrho, :null])
        @inbounds for i = 1:nznx
            mod.m[:vp][i] = χ(x[i], mod.ref[:vp], -1)
        end
        @inbounds for i = 1:nznx
            mod.m[:rho][i] = χ(x[nznx+i], mod.ref[:rho], -1)
        end
    elseif (fields == [:null, :χrho, :null])
        @inbounds for i = 1:nznx
            mod.m[:rho][i] = χ(x[i], mod.ref[:rho], -1)
        end
    elseif (fields == [:χvp, :null, :null])
        @inbounds for i = 1:nznx
            mod.m[:vp][i] = χ(x[i], mod.ref[:vp], -1)
        end
    else
        error("invalid fields")
    end
    return mod
end

"""
Input perturbations corresponding to `parameterization`.
Output perturbations of KI and rhoI, i.e., without the contrast.
`mod` is just used for reference values.
"""
function Base.copyto!(
    δxout::AbstractVector,
    δx::AbstractVector,
    mod::Medium,
    parameterization,
)
    nznx = prod(length.(mod.mgrid))
    fill!(δxout, 0.0)
    if (parameterization == [:χKI, :χrhoI, :null])
        @inbounds for i = 1:nznx
            δxout[i] = δx[i] * mod.ref[:KI]
            δxout[nznx+i] = δx[nznx+i] * mod.ref[:rhoI]
        end
    elseif (parameterization == [:χKI, :null, :null])
        @inbounds for i = 1:nznx
            δxout[i] = δx[i] * mod.ref[:KI]
        end
    elseif (parameterization == [:χvp, :χrho, :null])
        @inbounds for i = 1:nznx
            δxout[i] = (-1.0 * δx[nznx+i] - 2.0 * δx[i]) * mod.ref[:KI]
            δxout[nznx+i] = (-1.0 * δx[nznx+i]) * mod.ref[:rhoI]
        end
    elseif (parameterization == [:χvp, :null, :null])
        @inbounds for i = 1:nznx
            δxout[i] = -2.0 * δx[i] * mod.ref[:KI]
        end
    elseif (parameterization == [:null, :χrho, :null])
        @inbounds for i = 1:nznx
            δxout[i] = -1.0 * δx[i] * mod.ref[:rhoI]
        end
    end
end

#pert_of_KI(dvp, drho, vp0, rho0) = - 2. / (rho0 * vp0^3) * dvp - 1 / (rho0^2 * vp0^2) * drho
#pert_of_rhoI(drho, vp0, rho0) = -1. / (rho0^2) * drho

"""
Get other dependent model parameters of a seismic model
that are not present in `Medium`.

* `:rhoI` : inverse of density
* `:Zp` : P-wave impedance
* `:lambda` : first Lame parameter
* `:mu` : shear modulus
* `:K` : bulk modulus
* `:M` : P-wave modulus
"""
function Base.getindex(mod::Medium, s::Symbol)
    nznx = length(mod.m[:vp])
    # check if vp, vs, and rho are present
    @assert :vp ∈ names(mod.m)[1]
    @assert :rho ∈ names(mod.m)[1]
    vp = mod.m[:vp]
    rho = mod.m[:rho]
    if (:vs ∈ names(mod.m)[1])
        vs = mod.m[:vs]
    end
    if (s == :vp)
        return vp
    elseif (s == :rho)
        return rho
    elseif (s in [:vs, :χmu])
        @assert :vs ∈ names(mod.m)[1]
        if (s == :vs)
            return mod.m[:vs]
        end
    elseif (s == :Q)
        @assert :Q ∈ names(mod.m)[1]
        return mod.m[:Q]
    end

    # these are 3D arrays
    if (s in [:tau_epsilon, :tau_sigma])
        @assert :Q ∈ names(mod.m)[1] # need quality factor
        @assert :freqmin ∈ names(mod.fc)[1] # need minimum freq that is being modelled
        @assert :freqmax ∈ names(mod.fc)[1] # need max freq
        @assert :nsls ∈ names(mod.ic)[1] # need number of standard linear solids

        nz, nx = size(mod.m[:vp])
        # add derived fields, if not present
        (s ∉ names(mod.m3)[1]) &&
            (mod.m3 = vcat(mod.m3, NamedArray([zeros(mod.ic[:nsls], nz, nx)], [s])))
        x = mod.m3[s] # going to return this
    else
        # add derived fields, if not present
        (s ∉ names(mod.m)[1]) && (mod.m = vcat(mod.m, NamedArray([zero(vp)], [s])))
        x = mod.m[s] # going to return this
    end

    # output derived fields
    if (s == :rhoI)
        @inbounds for i = 1:nznx
            x[i] = inv(rho[i])
        end
    elseif (s == :lambda)
        if (:vs ∈ names(mod.m)[1])
            @inbounds for i = 1:nznx
                x[i] = (vp[i] * vp[i] - 2 * vs[i] * vs[i]) * rho[i]
            end
        else
            @inbounds for i = 1:nznx
                x[i] = vp[i] * vp[i] * rho[i]
            end
        end
    elseif (s == :M)
        @inbounds for i = 1:nznx
            x[i] = vp[i] * vp[i] * rho[i]
        end
    elseif (s == :χrhoI)
        @inbounds for i = 1:nznx
            x[i] = χ(inv(rho[i]), mod.ref[:rhoI], 1)
        end
    elseif (s == :χrho)
        @inbounds for i = 1:nznx
            x[i] = χ(rho[i], mod.ref[:rho], 1)
        end
    elseif (s == :χvp)
        @inbounds for i = 1:nznx
            x[i] = χ(vp[i], mod.ref[:vp], 1)
        end
    elseif (s == :χK)
        if (:vs ∈ names(mod.m)[1])
            @inbounds for i = 1:nznx
                x[i] =
                    χ((vp[i] * vp[i] - 4.0 / 3.0 * vs[i] * vs[i]) * rho[i], mod.ref[:K], 1)
            end
        else
            @inbounds for i = 1:nznx
                x[i] = χ(vp[i] * vp[i] * rho[i], mod.ref[:K], 1)
            end
        end
    elseif (s == :χM)
        @inbounds for i = 1:nznx
            x[i] = χ(vp[i] * vp[i] * rho[i], mod.ref[:K], 1)
        end
    elseif (s == :mu)
        @inbounds for i = 1:nznx
            x[i] = mod.m[:vs][i] * mod.m[:vs][i] * rho[i]
        end
    elseif (s == :χmu)
        @inbounds for i = 1:nznx
            x[i] = χ(mod.m[:vs][i] * mod.m[:vs][i] * rho[i], mod.ref[:mu], 1)
        end
    elseif (s == :χKI)
        if (:vs ∈ names(mod.m)[1])
            @inbounds for i = 1:nznx
                x[i] = χ(
                    inv((vp[i] * vp[i] - 4.0 / 3.0 * vs[i] * vs[i]) * rho[i]),
                    mod.ref[:KI],
                    1,
                )
            end
        else
            @inbounds for i = 1:nznx
                x[i] = χ(inv(vp[i] * vp[i] * rho[i]), mod.ref[:KI], 1)
            end
        end
    elseif (s == :χMI)
        @inbounds for i = 1:nznx
            x[i] = χ(inv(vp[i] * vp[i] * rho[i]), mod.ref[:KI], 1)
        end
    elseif (s == :KI)
        if (:vs ∈ names(mod.m)[1])
            @inbounds for i = 1:nznx
                x[i] = inv((vp[i] * vp[i] - 4.0 / 3.0 * vs[i] * vs[i]) * rho[i])
            end
        else
            @inbounds for i = 1:nznx
                x[i] = inv(vp[i] * vp[i] * rho[i])
            end
        end
    elseif (s == :K)
        if (:vs ∈ names(mod.m)[1])
            @inbounds for i = 1:nznx
                x[i] = (vp[i] * vp[i] - 4.0 / 3.0 * vs[i] * vs[i]) * rho[i]
            end
        else
            @inbounds for i = 1:nznx
                x[i] = vp[i] * vp[i] * rho[i]
            end
        end
    elseif (s == :Zp)
        @inbounds for i = 1:nznx
            x[i] = vp[i] * rho[i]
        end
    elseif (s == :tau_epsilon)
        # memory optimize later
        tau_sigma, tau_epsilon = get_relaxation_times(
            mod.m[:Q],
            mod.ic[:nsls],
            mod.fc[:freqmin],
            mod.fc[:freqmax],
        )
        copyto!(x, tau_epsilon)
    elseif (s == :tau_sigma)
        # memory optimize later
        tau_sigma, tau_epsilon = get_relaxation_times(
            mod.m[:Q],
            mod.ic[:nsls],
            mod.fc[:freqmin],
            mod.fc[:freqmax],
        )
        copyto!(x, tau_sigma)
    else
        error(string("invalid attrib\t", s))
    end
    return x
end

"""
Copy medium parameters from Medium `mod`, to vector `x`. Return `x`.
```julia
copyto!(x,mod,[:vp, :vs])
copyto!(x,mod,[:vp, :rho])
```
"""
function Base.copyto!(x::AbstractArray{T}, mod::Medium, fields::Vector{Symbol}) where {T}
    nznx = prod(length.(mod.mgrid))
    (length(x) ≠ (count(fields .≠ :null) * nznx)) && error("size x")
    i0 = 0
    for field in fields
        if (field ≠ :null)
            xx = view(x, i0+1:i0+nznx)
            m = Base.getindex(mod, field)
            copyto!(xx, m)
            i0 += nznx
        end
    end
    return x
end

"""
Return a vector of medium parameters `x`, from Medium `mod`.
```julia
x=vec(mod, field)
```
"""
function Base.vec(mod::Medium, field::Symbol)
    # allocate
    rout = zeros(length(mod.mgrid[1]), length(mod.mgrid[2]))
    copyto!(rout, mod, [field])
    return rout
end


text:
#using Images
#using Colors

struct Img
    path::String
    width::Int32
    height::Int32
    texture::Array{Cairo.ARGB32,1}
end

ARGB(r, g, b, a) = reinterpret(Cairo.ARGB32, UInt32(a)<<24 | UInt32(b)<<16 | UInt32(g)<<8 | UInt32(r))

function extractFromFile(path::String, filename::String)
    fh = open(joinpath(path, filename), "r")
    data = IOBuffer(read(fh))
    close(fh)
    width = read(data, Int32)
    height = read(data, Int32)
    hasAlpha = read(data, UInt8) == 1
    totalSize = width * height
    textureData = Array{Cairo.ARGB32,1}(totalSize)
    j = 1
    while j < totalSize
        repeat = read(data, UInt8)
        if hasAlpha
            alpha = read(data, UInt8)
            if alpha > 0
                textureData[j] = ARGB(read(data, UInt8), read(data, UInt8), read(data, UInt8), alpha)
            else
                textureData[j] = ARGB(UInt8(0), UInt8(0), UInt8(0), UInt8(0))
            end
        else
            textureData[j] = ARGB(read(data, UInt8), read(data, UInt8), read(data, UInt8), UInt8(0xFF))
        end
        if repeat > 1
            for i in 1:repeat - 1
                textureData[j + i] = deepcopy(textureData[j])
            end
        end
        j += repeat
    end
    Img(splitext(filename)[1], width, height, textureData)
end

function extractGraphics(celestedir::String)
    path = joinpath(celestedir, "Content", "Graphics", "Atlases")
    filepaths = String[]
    for (root, dirs, files) in walkdir(path), file in filter(f -> contains(f, ".data"), files)
        push!(filepaths, relpath(joinpath(root, file), path))
    end
    println("Extracting from $path - $(length(filepaths)) files")
    images = Img[]
    for file in filepaths
        println(file)
        push!(images, extractFromFile(path, file))
    end
    return images
end

function dumpSprites(path::String, destination::String=source)
    img = extractFromFile(dirname(path), basename(path))

    surface = CairoImageSurface(reshape(img.texture, (img.height, img.width)))
    open(io -> write_to_png(surface, io), destination, "w")
end

max_stars_repo_name: npdmfixup/Fuji

max_issues_repo_name: npdmfixup/Fuji

max_forks_repo_name: npdmfixup/Fuji


text:
export Vec
module Vec

# X[i] += a
function add_scal!{T}(X::Array{T}, a)
  leng = length(X)
  a = convert(eltype(X), a)
  @simd for i = 1:leng
    @inbounds X[i] += a
  end
end

# X[i] *= a
function mul_scal!{T}(X::Array{T}, a)
  leng = length(X)
  a = convert(eltype(X), a)
  @simd for i = 1:leng
    @inbounds X[i] *= a
  end
end

# X[i] *= Y[i]
function mul!{T}(X::Array{T}, Y::Array{T})
  leng = length(X)
  @simd for i = 1:leng
    @inbounds X[i] *= Y[i]
  end
end

# X[i] = X[i] / Y[i]
function div!{T}(X::Array{T}, Y::Array{T})
  leng = length(X)
  @simd for i = 1:leng
    @inbounds X[i] /= Y[i]
  end
end
# Y[i] = X[i] / Y[i]
function div2!{T}(X::Array{T}, Y::Array{T})
  leng = length(X)
  @simd for i = 1:leng
    @inbounds Y[i] = X[i] / Y[i]
  end
end

# X[i] = X[i]^p
function pow!{T}(X::Array{T}, p::Number)
  leng = length(X)
  @simd for i = 1:leng
    @inbounds X[i] = X[i]^p
  end
end

end # module Vec


text:
let
    ObjectID()
    RobotID()
    LocationID()
    ActionID()
    OperationID()

    get_id(ObjectID())

    @test ObjectID(1) < ObjectID(2)
    @test ObjectID(2) > ObjectID(1)
    @test ObjectID(1) + 1 == ObjectID(2)
    @test ObjectID(1) == ObjectID(2) - 1
end
let
    ROBOT_AT(RobotID(1),LocationID(2))
    ROBOT_AT(1,2)
    get_robot_id(ROBOT_AT(1,2))
end
let
    OBJECT_AT(1,2)
    OBJECT_AT(1,[2])
    OBJECT_AT(1,LocationID(2))
    OBJECT_AT(1,[LocationID(2)])

    get_object_id(OBJECT_AT(1,2))
    get_location_id(OBJECT_AT(1,2))
    get_location_ids(OBJECT_AT(1,2))
end
let
    op = Operation()
    set_precondition!(op,OBJECT_AT(1,2))
    @test TaskGraphs.get_dropoff(op,ObjectID(1)) == LocationID(2)
    get_location_id(op)
    get_operation_id(op)
    get_input_ids(op)
    get_output_ids(op)
    duration(op)
    preconditions(op)
    postconditions(op)
    get_id(op)
end
let
    GO()
    node = GO(1,2,3)
    get_initial_location_id(node)
    get_destination_location_id(node)
    get_robot_id(node)
    # GO does not have a single location, so this method should fail
    @test_throws MethodError get_location_id(node)
    @test_throws MethodError get_object_id(node)
end
let
    CARRY()
    node = CARRY(1,2,3,4)
    get_initial_location_id(node)
    get_destination_location_id(node)
    get_object_id(node)
    get_robot_id(node)
    # CARRY does not have a single location, so this method should fail
    @test_throws MethodError get_location_id(node)
end
let
    COLLECT()
    node = COLLECT(1,2,3)
    get_initial_location_id(node)
    get_destination_location_id(node)
    get_location_id(node)
    get_object_id(node)
    get_robot_id(node)
end
let
    DEPOSIT()
    node = DEPOSIT(1,2,3)
    get_initial_location_id(node)
    get_destination_location_id(node)
    get_location_id(node)
    get_object_id(node)
    get_robot_id(node)
end
let
    for node in [OBJECT_AT(1,2),ROBOT_AT(1,2),GO(),COLLECT(),CARRY(),DEPOSIT(),Operation()]
        required_predecessors(node)
        required_successors(node)
        eligible_predecessors(node)
        eligible_successors(node)
        resources_reserved(node)
        for node2 in [OBJECT_AT(1,2),ROBOT_AT(1,2),GO(),COLLECT(),CARRY(),DEPOSIT(),Operation()]
            matches_template(typeof(node),node2)
            matches_template((typeof(node),typeof(node2)),node2)
        end
    end
end
let
    for (T,nodes) in [
        AbstractPlanningPredicate=>[OBJECT_AT(1,2),ROBOT_AT(1,2),GO(),COLLECT(),CARRY(),DEPOSIT(),Operation()],
        AbstractSingleRobotAction=>[GO(),COLLECT(),CARRY(),DEPOSIT()],
        AbstractSingleRobotAction=>[CUB_GO(),CUB_COLLECT(),CUB_CARRY(),CUB_DEPOSIT()],
        ]
        for node in nodes
            @test matches_template(T,node)
            @test matches_template(T,typeof(node))
            @test matches_template((T,T),typeof(node))
            @test matches_template(node,node)
            @test matches_template((node,node),node)
            @test matches_template((node,GO()),node)
            @test matches_template(typeof(node),node)
            @test matches_template((typeof(node),typeof(node)),node)
            @test matches_template(T,ScheduleNode(RobotID(-1),node,PathSpec()))

            @test !matches_template(1,node)
            @test !matches_template((1,2),node)
            @test !matches_template(1,ScheduleNode(RobotID(-1),node,PathSpec()))
        end
    end
end
let
    @test split_node(GO(1,2,3),LocationID(4))[1] == GO(1,2,4)
    @test split_node(GO(1,2,3),LocationID(4))[2] == GO(1,4,3)
    @test split_node(CARRY(1,1,2,3),LocationID(4))[1] == CARRY(1,1,2,4)
    @test split_node(CARRY(1,1,2,3),LocationID(4))[2] == CARRY(1,1,4,3)
    @test split_node(COLLECT(1,2,3),LocationID(4))[1] == COLLECT(1,2,4)
    @test split_node(COLLECT(1,2,3),LocationID(4))[2] == COLLECT(1,2,4)
    @test split_node(DEPOSIT(1,2,3),LocationID(4))[1] == DEPOSIT(1,2,4)
    @test split_node(DEPOSIT(1,2,3),LocationID(4))[2] == DEPOSIT(1,2,4)
end
let
    for n in [GO(),CARRY(),COLLECT(),DEPOSIT()]
        @test TaskGraphs.replace_robot_id(n,RobotID(2)).r == RobotID(2)
    end
end
let
    n1 = GO(1,2,3)
    n2 = GO(-1,-1,4)
    n = align_with_predecessor(n2,n1)
    @test n.r == n1.r
    @test n.x1 == n1.x2
    @test n.x2 == n2.x2
end
let
    n1 = ROBOT_AT(1,2)
    n2 = GO(-1,-1,3)
    n = align_with_predecessor(n2,n1)
    @test n.r == n1.r
    @test n.x1 == n1.x
end
let
    n1 = DEPOSIT(1,2,3)
    n2 = GO(-1,3,4)
    n = align_with_predecessor(n2,n1)
    @test n.r == n1.r
    @test n.x1 == n1.x
end
let
    n1 = GO(1,2,4)
    n2 = COLLECT(-1,3,4)
    n = align_with_predecessor(n2,n1)
    @test n.r == n1.r
    @test n.x == n1.x2
end
let
    n1 = OBJECT_AT(3,4)
    n2 = COLLECT(-1,3,4)
    n = align_with_predecessor(n2,n1)
    @test n.o == n1.o
    @test n.x == n1.x[1]
end
let
    n1 = COLLECT(1,2,3)
    n2 = CARRY(-1,2,3,4)
    n = align_with_predecessor(n2,n1)
    @test n.o == n1.o
    @test n.x1 == n1.x
    @test n.r == n1.r
end
let
    n1 = CARRY(1,2,3,4)
    n2 = DEPOSIT(-1,2,4)
    n = align_with_predecessor(n2,n1)
    @test n.o == n1.o
    @test n.x == n1.x2
    @test n.r == n1.r
end
let
    TEAM_CARRY(instructions=[CARRY(1,2,3,4),CARRY(2,3,4,4)],shape=(2,1))

    sequence = [
        TEAM_COLLECT(instructions=[COLLECT(1,2,3),COLLECT(2,3,4)],shape=(2,1)),
        TEAM_CARRY(instructions=[CARRY(-1,2,3,4),CARRY(-1,3,4,5)],shape=(2,1)),
        TEAM_DEPOSIT(instructions=[DEPOSIT(-1,2,4),DEPOSIT(-1,3,5)],shape=(2,1)),
    ]
    @test all(map(p->get_id(get_robot_id(p)) == -1, sequence[2].instructions))
    sequence[2] = align_with_predecessor(sequence[2],sequence[1])
    @test all(map(p->get_id(get_robot_id(p)) != -1, sequence[2].instructions))

    @test all(map(p->get_id(get_robot_id(p)) == -1, sequence[3].instructions))
    sequence[3] = align_with_predecessor(sequence[3],sequence[2])
    @test all(map(p->get_id(get_robot_id(p)) != -1, sequence[3].instructions))
end
let
    @test validate_edge(ROBOT_AT(1,2),ROBOT_AT(1,2)) == false
    @test validate_edge(ROBOT_AT(1,2),GO(1,2,3))
    @test validate_edge(GO(1,2,3),GO(1,3,3))
    @test validate_edge(GO(1,2,3),COLLECT(1,4,3))
    @test validate_edge(COLLECT(1,2,3),CARRY(1,2,3,4))
    @test validate_edge(CARRY(1,2,3,4),COLLECT(1,2,3)) == false
    @test validate_edge(CARRY(1,2,3,4),DEPOSIT(1,2,4))
    @test validate_edge(DEPOSIT(1,2,4),CARRY(1,2,3,4)) == false
    @test validate_edge(DEPOSIT(1,2,4),GO(1,4,3))
    @test validate_edge(DEPOSIT(1,2,4),DEPOSIT(1,2,4))
    @test validate_edge(COLLECT(1,2,4),DEPOSIT(1,2,4))
    @test validate_edge(COLLECT(1,2,4),COLLECT(1,2,4))
end


text:
module JungleHelperDarkCrumbleBlock

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Crumble Blocks (Dark) (Jungle Helper)" => Ahorn.EntityPlacement(
        Maple.CrumbleBlock,
        "rectangle",
        Dict{String, Any}(
          "texture" => "JungleHelper/dark"
        )
    )
)

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
using DataFlow, DataFlow.Fuzz
using MacroTools, Lazy, Base.Test

import DataFlow: graphm, syntax, cse, dvertex, constant, prewalk

for nodes = 1:10, tries = 1:1_000

  dl = grow(DVertex, nodes)

  @test dl == @> dl syntax(bindconst = true) graphm

  @test copy(dl) == dl

  il = grow(IVertex, nodes)

  @test il == @> il DataFlow.dl() DataFlow.il()

  @test copy(il) == il == prewalk(identity, il)

end

@flow function recurrent(xs)
  hidden = σ( Wxh*xs + Whh*hidden + bh )
  σ( Wxy*x + Why*hidden + by )
end

@test @capture syntax(recurrent.output) begin
  h_Symbol = σ( Wxh*xs + Whh*h_Symbol + bh )
  σ( Wxy*x + Why*h_Symbol + by )
end

@flow function var(xs)
  mean = sum(xs)/length(xs)
  meansqr = sumabs2(xs)/length(xs)
  meansqr - mean^2
end

@test @capture syntax(var.output) begin
  sumabs2(xs)/length(xs) - (sum(xs) / length(xs)) ^ 2
end

@test contains(sprint(show, var),
               string(:(sumabs2(xs)/length(xs) - (sum(xs) / length(xs)) ^ 2)))

@test cse(var.output) == convert(IVertex, @flow begin
  n = length(xs)
  sumabs2(xs)/n - (sum(xs) / n) ^ 2
end)

let x = :(2+2)
  @test @flow(foo($x)) == dvertex(:foo, constant(x))
end


text:
import JLD2
using CUDA.CUDNN: cudnnConvolutionDescriptor

JLD2.writeas(::Type{cudnnConvolutionDescriptor}) = Nothing
JLD2.wconvert(::Type{Nothing}, c::cudnnConvolutionDescriptor) = nothing
JLD2.rconvert(::Type{cudnnConvolutionDescriptor}, j::Nothing) = cudnnConvolutionDescriptor(C_NULL) # TODO: fix this


text:
using LightGraphs
using MetaGraphs

function example_graph()
    fg = MetaDiGraph(SimpleDiGraph())
    nin = 3; nout = 2
    add_vertices!(fg, nin + nout + 3)
    types = [0*ones(nin); 2*ones(3); 1*ones(nout)]
    for i in eachindex(types)
        set_prop!(fg, i, :type, types[i])
    end
    set_prop!(fg, 4, :function, CGP.Config.f_sum)
    add_edge!(fg, Edge(1, 4))
    set_prop!(fg, 1, 4, :ci, 1)
    add_edge!(fg, Edge(2, 4))
    set_prop!(fg, 2, 4, :ci, 2)
    set_prop!(fg, 5, :function, CGP.Config.f_mult)
    add_edge!(fg, Edge(3, 5))
    set_prop!(fg, 3, 5, :ci, 1)
    add_edge!(fg, Edge(4, 5))
    set_prop!(fg, 4, 5, :ci, 2)
    set_prop!(fg, 6, :function, CGP.Config.f_aminus)
    add_edge!(fg, Edge(1, 6))
    set_prop!(fg, 1, 6, :ci, 1)
    add_edge!(fg, Edge(3, 6))
    set_prop!(fg, 3, 6, :ci, 2)
    add_edge!(fg, Edge(5, 7))
    set_prop!(fg, 5, 7, :ci, 0)
    add_edge!(fg, Edge(6, 8))
    set_prop!(fg, 6, 8, :ci, 0)
    fg
end


text:
# Functions to compute transformations between known datums

#
# ITRF -> ITRF transforms
#

"""
    ITRF_from_ITRF(ITRF_output_year, ITRF_input_year, epoch)

Return a `Transformation` converting ECEF points from one ITRF realization to another to GDA94.  Datum
shift parameters are taken from [1], supporting `ITRF_year` 2014, 2008, 2005, 2000, 1997, 1996
1997, 1996, 1994, 1993, 1992, 1991, 1990, 1989, 1988.  `epoch` is the `Date` (or `DateTime`) of interest at which the
input `ECEF` coordinates were measured in ITRF.

[1] http://itrf.ign.fr/trans_para.php

TODO: We don't yet support `epoch` varying per input point, but there should be
a `Transformation` object for this at some stage.

"""

function ITRF_from_ITRF(ITRF_output_year, ITRF_input_year, epoch)

    # Values comes from http://itrf.ign.fr/trans_para.php
    #
    # (Unit key: ppb = parts per billion; mas = milli arc seconds; yr = year)
    #
    #          Values at 2000 reference epoch                                 Derivatives                                          ReferenceEpoch
    #  Name    tx      ty      tz     D       rx       ry      rz             tx      ty    tz     D      rx      ry      rz
    #  Units   mm      mm      mm     ppb     mas      mas     mas            mm/yr  mm/yr  mm/yr  ppb/yr mas/yr  mas/yr  mas/yr
    #
    # dict key is (input, output)
    table = Dict(
        # 2014 ->
        (2014, 2005) => ([ 2.6000,  1.0000, -2.3000,  0.9200,  0.0000,  0.0000,  0.0000], [ 0.3000,  0.0000, -0.1000,  0.0300,  0.0000,  0.0000,  0.0000], 2010),
        (2014, 2000) => ([ 0.7000,  1.2000, -26.1000,  2.1200,  0.0000,  0.0000,  0.0000], [ 0.1000,  0.1000, -1.9000,  0.1100,  0.0000,  0.0000,  0.0000], 2010),
        (2014, 1997) => ([ 7.4000, -0.5000, -62.8000,  3.8000,  0.0000,  0.0000,  0.2600], [ 0.1000, -0.5000, -3.3000,  0.1200,  0.0000,  0.0000,  0.0200], 2010),
        (2014, 1996) => ([ 7.4000, -0.5000, -62.8000,  3.8000,  0.0000,  0.0000,  0.2600], [ 0.1000, -0.5000, -3.3000,  0.1200,  0.0000,  0.0000,  0.0200], 2010),
        (2014, 1994) => ([ 7.4000, -0.5000, -62.8000,  3.8000,  0.0000,  0.0000,  0.2600], [ 0.1000, -0.5000, -3.3000,  0.1200,  0.0000,  0.0000,  0.0200], 2010),
        (2014, 1993) => ([-50.4000,  3.3000, -60.2000,  4.2900, -2.8100, -3.3800,  0.4000], [-2.8000, -0.1000, -2.5000,  0.1200, -0.1100, -0.1900,  0.0700], 2010),
        (2014, 1992) => ([ 15.4000,  1.5000, -70.8000,  3.0900,  0.0000,  0.0000,  0.2600], [ 0.1000, -0.5000, -3.3000,  0.1200,  0.0000,  0.0000,  0.0200], 2010),
        (2014, 1991) => ([ 27.4000,  15.5000, -76.8000,  4.4900,  0.0000,  0.0000,  0.2600], [ 0.1000, -0.5000, -3.3000,  0.1200,  0.0000,  0.0000,  0.0200], 2010),
        (2014, 1990) => ([ 25.4000,  11.5000, -92.8000,  4.7900,  0.0000,  0.0000,  0.2600], [ 0.1000, -0.5000, -3.3000,  0.1200,  0.0000,  0.0000,  0.0200], 2010),
        (2014, 1989) => ([ 30.4000,  35.5000, -130.8000,  8.1900,  0.0000,  0.0000,  0.2600], [ 0.1000, -0.5000, -3.3000,  0.1200,  0.0000,  0.0000,  0.0200], 2010),
        (2014, 1988) => ([ 25.4000, -0.5000, -154.8000,  11.2900,  0.1000,  0.0000,  0.2600], [ 0.1000, -0.5000, -3.3000,  0.1200,  0.0000,  0.0000,  0.0200], 2010),
        # 2008 ->
        (2008, 2005) => ([-2.0000, -0.9000, -4.7000,  0.9400,  0.0000,  0.0000,  0.0000], [ 0.3000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000], 2000),
        (2008, 2000) => ([-1.9000, -1.7000, -10.5000,  1.3400,  0.0000,  0.0000,  0.0000], [ 0.1000,  0.1000, -1.8000,  0.0800,  0.0000,  0.0000,  0.0000], 2000),
        (2008, 1997) => ([ 4.8000,  2.6000, -33.2000,  2.9200,  0.0000,  0.0000,  0.0600], [ 0.1000, -0.5000, -3.2000,  0.0900,  0.0000,  0.0000,  0.0200], 2000),
        (2008, 1996) => ([ 4.8000,  2.6000, -33.2000,  2.9200,  0.0000,  0.0000,  0.0600], [ 0.1000, -0.5000, -3.2000,  0.0900,  0.0000,  0.0000,  0.0200], 2000),
        (2008, 1994) => ([ 4.8000,  2.6000, -33.2000,  2.9200,  0.0000,  0.0000,  0.0600], [ 0.1000, -0.5000, -3.2000,  0.0900,  0.0000,  0.0000,  0.0200], 2000),
        (2008, 1993) => ([-24.0000,  2.4000, -38.6000,  3.4100, -1.7100, -1.4800, -0.3000], [-2.8000, -0.1000, -2.4000,  0.0900, -0.1100, -0.1900,  0.0700], 2000),
        (2008, 1992) => ([ 12.8000,  4.6000, -41.2000,  2.2100,  0.0000,  0.0000,  0.0600], [ 0.1000, -0.5000, -3.2000,  0.0900,  0.0000,  0.0000,  0.0200], 2000),
        (2008, 1991) => ([ 24.8000,  18.6000, -47.2000,  3.6100,  0.0000,  0.0000,  0.0600], [ 0.1000, -0.5000, -3.2000,  0.0900,  0.0000,  0.0000,  0.0200], 2000),
        (2008, 1990) => ([ 22.8000,  14.6000, -63.2000,  3.9100,  0.0000,  0.0000,  0.0600], [ 0.1000, -0.5000, -3.2000,  0.0900,  0.0000,  0.0000,  0.0200], 2000),
        (2008, 1989) => ([ 27.8000,  38.6000, -101.2000,  7.3100,  0.0000,  0.0000,  0.0600], [ 0.1000, -0.5000, -3.2000,  0.0900,  0.0000,  0.0000,  0.0200], 2000),
        (2008, 1988) => ([ 22.8000,  2.6000, -125.2000,  10.4100,  0.1000,  0.0000,  0.0600], [ 0.1000, -0.5000, -3.2000,  0.0900,  0.0000,  0.0000,  0.0200], 2000),
        # 2005 ->
        (2005, 2000) => ([0.1 -0.8 -5.8 0.40 0.000 0.000 0.000], [-0.2 0.1 -1.8 0.08 0.000 0.000 0.000], 2000),
        # 2000 ->
        (2000, 1997) => ([ 0.6700,  0.6100, -1.8500,  1.5500,  0.0000,  0.0000,  0.0000], [ 0.0000, -0.0600, -0.1400,  0.0100,  0.0000,  0.0000,  0.0200], 1997),
        (2000, 1996) => ([ 0.6700,  0.6100, -1.8500,  1.5500,  0.0000,  0.0000,  0.0000], [ 0.0000, -0.0600, -0.1400,  0.0100,  0.0000,  0.0000,  0.0200], 1997),
        (2000, 1994) => ([ 0.6700,  0.6100, -1.8500,  1.5500,  0.0000,  0.0000,  0.0000], [ 0.0000, -0.0600, -0.1400,  0.0100,  0.0000,  0.0000,  0.0200], 1997),
        (2000, 1993) => ([ 1.2700,  0.6500, -2.0900,  1.9500, -0.3900,  0.8000, -1.1400], [-0.2900, -0.0200, -0.0600,  0.0100, -0.1100, -0.1900,  0.0700], 1988),
        (2000, 1992) => ([ 1.4700,  1.3500, -1.3900,  0.7500,  0.0000,  0.0000, -0.1800], [ 0.0000, -0.0600, -0.1400,  0.0100,  0.0000,  0.0000,  0.0200], 1988),
        (2000, 1991) => ([ 2.6700,  2.7500, -1.9900,  2.1500,  0.0000,  0.0000, -0.1800], [ 0.0000, -0.0600, -0.1400,  0.0100,  0.0000,  0.0000,  0.0200], 1988),
        (2000, 1990) => ([ 2.4700,  2.3500, -3.5900,  2.4500,  0.0000,  0.0000, -0.1800], [ 0.0000, -0.0600, -0.1400,  0.0100,  0.0000,  0.0000,  0.0200], 1988),
        (2000, 1989) => ([ 2.9700,  4.7500, -7.3900,  5.8500,  0.0000,  0.0000, -0.1800], [ 0.0000, -0.0600, -0.1400,  0.0100,  0.0000,  0.0000,  0.0200], 1988),
        (2000, 1988) => ([ 2.4700,  1.1500, -9.7900,  8.9500,  0.1000,  0.0000, -0.1800], [ 0.0000, -0.0600, -0.1400,  0.0100,  0.0000,  0.0000,  0.0200], 1988)
    )
    key = (max(ITRF_output_year, ITRF_input_year), min(ITRF_output_year, ITRF_input_year))
    if !haskey(table, key)
        throw(ErrorException("No ITRF$(ITRF_input_year) to ITRF$(ITRF_output_year) transformation available"))
    end
    params, rates, reference_epoch = table[key]
    reference_epoch = DateTime(reference_epoch,1,1)

    # Fractional years since reference epoch.  Seems a bit unnecessarily cumbersome!
    millisecs_per_year = 365.25 * Dates.value(Dates.Millisecond(Dates.Day(1)))
    dt = Dates.value(Dates.Millisecond(convert(DateTime, epoch) - reference_epoch)) / millisecs_per_year

    # Convert units to meters and radians
    mas2rad = deg2rad(1e-3/(60*60))
    unitconv = [1e-3, 1e-3, 1e-3, 1e-9, mas2rad, mas2rad, mas2rad]
    Tx,Ty,Tz, D, Rx,Ry,Rz = unitconv .* (params .+ rates*dt)

    # assemble the backwards in time version
    M = @fsa [1+D  Rz  -Ry;
             -Rz   1+D  Rx;
              Ry  -Rx   1+D]
    T = Vec(Tx,Ty,Tz)

    # want the inverse to go forward in time
    if (ITRF_output_year > ITRF_input_year)
        M = inv(M)
        T = -M * T
    end
    return (M, T)
end

#
# GDA94 related transforms
#

"""
    GDA94_from_ITRF(ITRF_year, epoch)

Return a `Transformation` converting ECEF points from ITRF to GDA94.  Datum
shift parameters are taken from [1], supporting `ITRF_year` 2008, 2005, 2000,
1997, 1996.  `epoch` is the `Date` (or `DateTime`) of interest at which the
input `ECEF` coordinates were measured in ITRF.

[1] J. Dawson and A. Woods, "ITRF to GDA94 coordinate transforms",
    Journal of Applied Geodesy, 4, p. 189 (2010).

TODO: We don't yet support `epoch` varying per input point, but there should be
a `Transformation` object for this at some stage.
"""
function GDA94_from_ITRF(ITRF_year, epoch)
    # ITRF transformation parameters from
    # J. Dawson and A. Woods, "ITRF to GDA94 coordinate transforms",
    # Journal of Applied Geodesy, 4, p. 189 (2010):
    #
    # These are intentionally kept exactly the same as the paper for
    # consistency and ease of checking.
    #
    # (Unit key: ppb = parts per billion; mas = milli arc seconds; yr = year)
    #
    #          Values at 1994 reference epoch                                 Derivatives
    #  Name    tx      ty      tz     sc      rx       ry      rz             tx      ty    tz     sc     rx      ry      rz
    #  Units   mm      mm      mm     ppb     mas      mas     mas            mm/yr  mm/yr  mm/yr  ppb/yr mas/yr  mas/yr  mas/yr
    #
    table = Dict(
        2008 => ([-84.68, -19.42,  32.01, 9.710, -0.4254,  2.2578, 2.4015], [  1.42,  1.34,  0.90, 0.109, 1.5461, 1.1820, 1.1551]),
        2005 => ([-79.73,  -6.86,  38.03, 6.636, -0.0351,  2.1211, 2.1411], [  2.25, -0.62, -0.56, 0.294, 1.4707, 1.1443, 1.1701]),
        2000 => ([-45.91, -29.85, -20.37, 7.070, -1.6705,  0.4594, 1.9356], [ -4.66,  3.55, 11.24, 0.249, 1.7454, 1.4868, 1.2240]),
        1997 => ([-14.63, -27.62, -25.32, 6.695, -1.7893, -0.6047, 0.9962], [ -8.60,  0.36, 11.25, 0.007, 1.6394, 1.5198, 1.3801]),
        1996 => ([ 24.54, -36.43, -68.12, 6.901, -2.7359, -2.0431, 0.3731], [-21.80,  4.71, 26.27, 0.388, 2.0203, 2.1735, 1.6290]),
    )
    reference_epoch = DateTime(1994,1,1)
    (2008, 1994) => ([ 4.8000,  2.6000, -33.2000,  2.9200,  0.0000,  0.0000,  0.0600], [ 0.1000, -0.5000, -3.2000,  0.0900,  0.0000,  0.0000,  0.0200], 2000),

    if !haskey(table, ITRF_year)
        throw(ErrorException("No ITRF to GDA94 parameters for ITRF year $ITRF_year"))
    end
    params, rates = table[ITRF_year]

    # Fractional years since reference epoch.  Seems a bit unnecessarily cumbersome!
    millisecs_per_year = 365.25 * Dates.value(Dates.Millisecond(Dates.Day(1)))
    dt = Dates.value(Dates.Millisecond(convert(DateTime, epoch) - reference_epoch)) / millisecs_per_year

    # Convert units to meters and radians
    mas2rad = deg2rad(1e-3/(60*60))
    unitconv = [1e-3, 1e-3, 1e-3, 1e-9, mas2rad, mas2rad, mas2rad]
    Tx,Ty,Tz, Sc, Rx,Ry,Rz = unitconv .* (params .+ rates*dt)

    M = (1+Sc) * @fsa [1.0  Rz  -Ry;
                      -Rz  1.0   Rx;
                       Ry  -Rx  1.0]

    (M, Vec(Tx,Ty,Tz))
end

"""
    ITRF_from_GDA94(ITRF_year, epoch)

The inverse of `GDA94_from_ITRF()`: Return a `Transformation` converting ECEF
points from ITRF to GDA94.
"""
function ITRF_from_GDA94(ITRF_year, epoch)
    M, T = (GDA94_from_ITRF(ITRF_year, epoch))
    Mi = inv(M)
    return Mi, -Mi * T
end

#
# ETRS related transforms
#

"""
    ETRF_from_ITRF(ITRF_year, epoch)

Return a `Transformation` converting ECEF points from ITRF to ETRF2000.  Datum
shift parameters are taken from [1], supporting `ITRF_year` 2008, 2005, 2000,
1997, 1996, 1994, 1993, 1992, 1991, 1990, 1989.  `epoch` is the `Date` (or `DateTime`) of interest at which the
input `ECEF` coordinates were measured in ITRF.

[1] C. Boucher and Z. Altamimi, "Memo : Specifications for reference frame fixing in the analysis of a
    EUREF GPS campaign" Version 8

"""
function ETRF_from_ITRF(ETRF_year, ITRF_year, epoch)
    (M1, T1) = ITRF_from_ITRF(ETRF_year, ITRF_year, epoch)
    (M2, T2) = ETRF_from_ITRF(ETRF_year, epoch)
    (M2 * M1, M2 * T1 + T2) # compose
end

function ETRF_from_ITRF(Year, epoch)

    # ITRF transformation parameters from
    # C. Boucher and Z. Altamimi, "Memo : Specifications for reference frame fixing in the analysis of a
    # EUREF GPS campaign" Version 8
    #
    # These are intentionally kept exactly the same as the paper for
    # consistency and ease of checking.
    #
    # (Unit key: ppb = parts per billion; mas = milli arc seconds; yr = year)
    #
    #          Values at 1989 reference epoch                                 Derivatives
    #  Name    tx      ty      tz     D       rx       ry      rz             tx      ty    tz     D      rx      ry      rz
    #  Units   cm      cm      cm     1e-8    mas      mas     mas            cm/yr   cm/yr cm/yr  1e-8/yr mas/yr  mas/yr  mas/yr
    #
    table = Dict(
        2005 => ([ 5.6000,  4.8000, -3.7000,  0.0000,  0.0000,  0.0000,  0.0000], [ 0.0000,  0.0000,  0.0000,  0.0000,  0.0540,  0.5180, -0.7810]),
        2000 => ([ 5.4000,  5.1000, -4.8000,  0.0000,  0.0000,  0.0000,  0.0000], [ 0.0000,  0.0000,  0.0000,  0.0000,  0.0810,  0.4900, -0.7920]),
        1997 => ([ 4.1000,  4.1000, -4.9000,  0.0000,  0.0000,  0.0000,  0.0000], [ 0.0000,  0.0000,  0.0000,  0.0000,  0.2000,  0.5000, -0.6500]),
        1996 => ([ 4.1000,  4.1000, -4.9000,  0.0000,  0.0000,  0.0000,  0.0000], [ 0.0000,  0.0000,  0.0000,  0.0000,  0.2000,  0.5000, -0.6500]),
        1994 => ([ 4.1000,  4.1000, -4.9000,  0.0000,  0.0000,  0.0000,  0.0000], [ 0.0000,  0.0000,  0.0000,  0.0000,  0.2000,  0.5000, -0.6500]),
        1993 => ([ 1.9000,  5.3000, -2.1000,  0.0000,  0.0000,  0.0000,  0.0000], [ 0.0000,  0.0000,  0.0000,  0.0000,  0.3200,  0.7800, -0.6700]),
        1992 => ([ 3.8000,  4.0000, -3.7000,  0.0000,  0.0000,  0.0000,  0.0000], [ 0.0000,  0.0000,  0.0000,  0.0000,  0.2100,  0.5200, -0.6800]),
        1991 => ([ 2.1000,  2.5000, -3.7000,  0.0000,  0.0000,  0.0000,  0.0000], [ 0.0000,  0.0000,  0.0000,  0.0000,  0.2100,  0.5200, -0.6800]),
        1990 => ([ 1.9000,  2.8000, -2.3000,  0.0000,  0.0000,  0.0000,  0.0000], [ 0.0000,  0.0000,  0.0000,  0.0000,  0.1100,  0.5700, -0.7100]),
        1989 => ([ 0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000,  0.0000], [ 0.0000,  0.0000,  0.0000,  0.0000,  0.1100,  0.5700, -0.7100])
    )
    reference_epoch = DateTime(1989,1,1)

    if !haskey(table, Year)
        throw(ErrorException("No $(src) to $(dest) parameters available"))
    end
    params, rates = table[Year]

    # Fractional years since reference epoch.  Seems a bit unnecessarily cumbersome!
    millisecs_per_year = 365.25 * Dates.value(Dates.Millisecond(Dates.Day(1)))
    dt = Dates.value(Dates.Millisecond(convert(DateTime, epoch) - reference_epoch)) / millisecs_per_year

    # Convert units to meters and radians
    mas2rad = deg2rad(1e-3/(60*60))
    unitconv = [1e-2, 1e-2, 1e-2, 1e-8, mas2rad, mas2rad, mas2rad]
    Tx,Ty,Tz, D, Rx,Ry,Rz = unitconv .* (params .+ rates*dt)
    M =  @fsa [1+D  Rz  -Ry;
              -Rz   1+D  Rx;
               Ry  -Rx   1+D]
    (M, Vec(Tx,Ty,Tz))
end

function ETRF2000_from_ITRF(ITRF_year, epoch)
    # ITRF transformation parameters from
    # C. Boucher and Z. Altamimi, "Memo : Specifications for reference frame fixing in the analysis of a
    # EUREF GPS campaign" Version 8
    #
    # These are intentionally kept exactly the same as the paper for
    # consistency and ease of checking.
    #
    # (Unit key: ppb = parts per billion; mas = milli arc seconds; yr = year)
    #
    #          Values at 2000 reference epoch                                 Derivatives
    #  Name    tx      ty      tz     D       rx       ry      rz             tx      ty    tz     D      rx      ry      rz
    #  Units   mm      mm      mm     ppb     mas      mas     mas            mm/yr  mm/yr  mm/yr  ppb/yr mas/yr  mas/yr  mas/yr
    #
    table = Dict(
        2008 => ([52.1, 49.3, -58.5,  1.34, 0.891, 5.390, -8.712], [0.1,  0.1, -1.8,  0.08, 0.081, 0.490, -0.792]),
        2005 => ([54.1, 50.2, -53.8,  0.40, 0.891, 5.390, -8.712], [-0.2, 0.1, -1.8,  0.08, 0.081, 0.490, -0.792]),
        2000 => ([54.0, 51.0, -48.0,  0.00, 0.891, 5.390, -8.712], [0.0,  0.0,  0.0,  0.00, 0.081, 0.490, -0.792]),
        1997 => ([47.3, 46.7, -25.3, -1.58, 0.891, 5.390, -8.772], [0.0,  0.6,  1.4, -0.01, 0.081, 0.490, -0.812]),
        1996 => ([47.3, 46.7, -25.3, -1.58, 0.891, 5.390, -8.772], [0.0,  0.6,  1.4, -0.01, 0.081, 0.490, -0.812]),
        1994 => ([47.3, 46.7, -25.3, -1.58, 0.891, 5.390, -8.772], [0.0,  0.6,  1.4, -0.01, 0.081, 0.490, -0.812]),
        1993 => ([76.1, 46.9, -19.9, -2.07, 2.601, 6.870, -8.412], [2.9,  0.2,  0.6, -0.01, 0.191, 0.680, -0.862]),
        1992 => ([39.3, 44.7, -17.3, -0.87, 0.891, 5.390, -8.772], [0.0,  0.6,  1.4, -0.01, 0.081, 0.490, -0.812]),
        1991 => ([27.3, 30.7, -11.3, -2.27, 0.891, 5.390, -8.772], [0.0,  0.6,  1.4, -0.01, 0.081, 0.490, -0.812]),
        1990 => ([29.3, 34.7,  4.7,  -2.57, 0.891, 5.390, -8.772], [0.0,  0.6,  1.4, -0.01, 0.081, 0.490, -0.812]),
        1989 => ([24.3, 10.7,  42.7, -5.97, 0.891, 5.390, -8.772], [0.0,  0.6,  1.4, -0.01, 0.081, 0.490, -0.812])
    )
    reference_epoch = DateTime(2000,1,1)

    if !haskey(table, ITRF_year)
        throw(ErrorException("No ITRF to ETRF2000 parameters for ITRF year $ITRF_year"))
    end
    params, rates = table[ITRF_year]

    # Fractional years since reference epoch.  Seems a bit unnecessarily cumbersome!
    millisecs_per_year = 365.25 * Dates.value(Dates.Millisecond(Dates.Day(1)))
    dt = Dates.value(Dates.Millisecond(convert(DateTime, epoch) - reference_epoch)) / millisecs_per_year

    # Convert units to meters and radians
    mas2rad = deg2rad(1e-3/(60*60))
    unitconv = [1e-3, 1e-3, 1e-3, 1e-9, mas2rad, mas2rad, mas2rad]
    Tx,Ty,Tz, D, Rx,Ry,Rz = unitconv .* (params .+ rates*dt)

    M = @fsa [1+D  Rz  -Ry;
             -Rz   1+D  Rx;
              Ry  -Rx   1+D]

    (M, Vec(Tx,Ty,Tz))
end

"""
    ITRF_from_ETRF2000(ITRF_year, epoch)

The inverse of `ETRF2000_from_ITRF()`: Return a `Transformation` converting ECEF
points from ITRF to GDA94.
"""
function ITRF_from_ETRF(ITRF_year, ETRF_year, epoch)
    M, T = (ETRF_from_ITRF(ETRF_year, ITRF_year, epoch))
    Mi = inv(M)
    return Mi, -Mi * T
end

function ITRF_from_ETRF2000(ITRF_year, epoch)
    M, T = (ETRF2000_from_ITRF(ITRF_year, epoch))
    Mi = inv(M)
    return Mi, -Mi * T
end

#
# Helper functions
#

# parse the itrf conversion files and return a dict
# example http://itrf.ign.fr/doc_ITRF/Transfo-ITRF2014_ITRFs.txt
# hope they keep using the same format
function parse_itrf_src(file)

    # (year_out, year_in) => (params, rates, epeoch)
    dict = Dict{NTuple{2,Int}, Tuple{Vector{Float64}, Vector{Float64}, Int}}()

    # read it all
    local lines
    open(file) do fid

        # grab the realization year
        header = readline(fid)
        src_year = match(r"(?<=ITRF)\d+", header).match
        src_year = (length(src_year) < 4 ? 1900 : 0) + parse(Int, src_year)

        # now read till we get into the parameters part
        line = chomp(readline(fid))
        done = false
        while !ismatch(r"(?<=ITRF)\d+", line)
            line = chomp(readline(fid))
        end

        # and loop till thge end of the table
        while ismatch(r"(?<=ITRF)\d+", line)

            # read the year
            dest_year = match(r"(?<=ITRF)\d+", line).match
            dest_year = (length(dest_year) < 4 ? 1900 : 0) + parse(Int, dest_year)

            # read params
            params = Vector{Float64}(0)
            for str in eachmatch(r"(?<=\s)[-+]*\d+\.\d+", line); push!(params, parse(Float64, str.match)); end
            epoch = Int(pop!(params))
            @assert length(params) == 7

            # read params
            line = chomp(readline(fid))
            rates = Vector{Float64}(0)
            for str in eachmatch(r"(?<=\s)[-+]*\d+\.\d+", line); push!(rates, parse(Float64, str.match)); end
            @assert length(rates) == 7

            # add it to the dict
            dict[(src_year, dest_year)] = (params, rates, epoch)
            @printf("(%i, %i) => ([% 0.4f, % 0.4f, % 0.4f, % 0.4f, % 0.4f, % 0.4f, % 0.4f], [% 0.4f, % 0.4f, % 0.4f, % 0.4f, % 0.4f, % 0.4f, % 0.4f], %i),\n",
                     src_year, dest_year, params..., rates..., epoch)

            # move along
            line = chomp(readline(fid))
        end
    end
    return dict
end


text:
# Autogenerated wrapper script for Gloo_jll for x86_64-linux-gnu-cxx11
export libgloo

JLLWrappers.@generate_wrapper_header("Gloo")
JLLWrappers.@declare_library_product(libgloo, "libgloo.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libgloo,
        "lib/libgloo.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using DiscreteFunctions
using Test

@testset "Segment"
    s1 = Segment(3)
    s2 = Segment(1, 3)
    @test iterator(s1) = 0 : 2
    @test iterator(s2) = 1 : 3
    @test first(s1) == 0
    @test last(s2) == 3
    @test length(s1) == 3
    @test length(s2) == 3
    # iterators
    s = 0
    for x in s1
        s += x
    end
    @test s == 2
    s = 0
    for x in s2
        s += x
    end
    @test s == 6
    @test ndims(s1) == 1
    @test size(s1) == (3,)
end

@testset "BooleanCube"
    bc = BooleanCube(5)
    @test length(bc) == 2^5
    # iterators
    s = 0
    for x in bc
        s += value(x)
    end
    @test s == 496
end

@testset "Direct product"
    bc = BooleanCube(5)
    s1 = Segment(3)
    dp = DirectProduct(bc, s1)
    @test length(bc) == 2^5 * 3
    s = 0
    for x in bc
        s += value(x)
    end
    @test s == 496
end

max_stars_repo_name: kirtsar/DiscreteFunctions

max_issues_repo_name: kirtsar/DiscreteFunctions


text:
# make a configuration broadcastable
Broadcast.broadcastable( sm::SurrogateModel ) = Ref(sm);
Broadcast.broadcastable( sc::SurrogateConfig ) = Ref(sc);

# Methods to be implemented by each type inheriting from SurrogateConfig
max_evals( :: SurrogateConfig ) = typemax(Int) :: Int;

# return data that is stored in iter data in each iteration
# get_saveable(::SurrogateMeta) = nothing :: Union{Nothing, <:SurrogateMeta};

fully_linear( :: SurrogateModel ) = false :: Bool;

# can objective functions with same configuration types be combined
# to a new vector objective?
combinable( :: SurrogateConfig ) = false :: Bool

# TODO make combinable bi-variate to check for to concrete configs if they are combinable

## TODO: make `prepare_init_model` and `_init_model` have a `ensure_fully_linear` kwarg too
function prepare_init_model( ::SurrogateConfig, :: AbstractObjective, :: AbstractMOP,
    :: AbstractIterData, ::AbstractDB, :: AbstractConfig; kwargs... ) :: SurrogateMeta
    nothing
end

function _init_model( ::SurrogateConfig, :: AbstractObjective, :: AbstractMOP,
    :: AbstractIterData, ::AbstractDB, :: AbstractConfig, :: SurrogateMeta; kwargs... ) :: Tuple{<:SurrogateModel,<:SurrogateMeta}
    nothing
end

## TODO: Allow to pass a SurrogateConfig here as well. (ATM use `model_cfg(objf)`) #src
## In general, the function signatures are somewhat messy. We should unify them a bit. #src
function update_model( :: SurrogateModel, :: AbstractObjective, :: SurrogateMeta, :: AbstractMOP,
    :: AbstractIterData, :: AbstractDB, :: AbstractConfig; kwargs... ) :: Tuple{<:SurrogateModel,<:SurrogateMeta}
    nothing
end

eval_models( :: SurrogateModel, ::Vec ) ::Vec = nothing
get_gradient( :: SurrogateModel, ::Vec, :: Int ) :: Vec = nothing
get_jacobian( :: SurrogateModel, :: Vec ) :: Mat = nothing

# DEFAULTS

get_saveable_type( :: SurrogateMeta ) = Nothing
get_saveable( :: SurrogateMeta ) = nothing

prepare_update_model( mod, objf, meta, mop, iter_data, db, algo_config; kwargs...) = meta
prepare_improve_model( mod, objf, meta, mop, iter_data, db, algo_config; kwargs...) = meta

# overwrite if possible, this is inefficient:
eval_models( sm :: SurrogateModel, x̂ :: Vec, ℓ :: Int) = eval_models(sm, x̂)[ℓ]

function improve_model( mod :: SurrogateModel, objf:: AbstractObjective, meta :: SurrogateMeta,
    mop :: AbstractMOP, id :: AbstractIterData, db :: AbstractDB, ac :: AbstractConfig;
    ensure_fully_linear = false)
    return mod, meta
end

# check if surrogate configurations are equal (only really needed if combinable)
function Base.:(==)( cfg1 :: T, cfg2 :: T ) where T <: SurrogateConfig
    all( getfield(cfg1, fname) == getfield(cfg2, fname) for fname ∈ fieldnames(T) )
end

function Base.:(==)( cfg1 :: T, cfg2 :: F ) where {T <: SurrogateConfig, F<:SurrogateConfig}
    false
end

# only needed if combinable
combine( cfg :: SurrogateConfig, :: SurrogateConfig  ) = cfg

## derived

@doc """
Return a function handle to be used with `NLopt` for output `ℓ` of `model`.
That is, if `model` is a surrogate for two scalar objectives, then `ℓ` must
be either 1 or 2.
"""
function _get_optim_handle( model :: SurrogateModel, ℓ :: Int )
    # Return an anonymous function that modifies the gradient if present
    function (x :: Vec, g :: Vec)
        if !isempty(g)
            g[:] = get_gradient( model, x, ℓ)
        end
        return eval_models( model, x, ℓ)
    end
end


text:
import Logging: Info,
                default_logcolor,
                min_enabled_level,
                shouldlog,
                handle_message,
                termlength,
                showvalue

import Base.CoreLogging: logmsg_code,
                @_sourceinfo,
                _min_enabled_level,
                current_logger_for_env,
                logging_error
import Base.string_with_env
const Wslog = LogLevel(50)
"""
Differences to stdlib/Logging/ConsoleLogger:

    - default timestamp on logging messages (except @info)
    - a 'shouldlog' function can be passed in. The `shouldlog_default` function filters
        on HTTP.Servers messages as well as on message_limits
    - :wslog => true flag which may be used for context-sensitive output
        from 'show' methods. This means a user can define 'show' methods
        which are used with this logger without affecting the behaviour
        defined in other modules.
    - :limited => true is included in the default IOContext. Keyword: show_limited
    - string_with_env_ws is exported for easy overloading on specific types
    - @info, @debug, @warn etc. will splat the first argument if it's a tuple arguments, e.g.

    julia> var = "a"
    "a"
    julia> @info (1, var)
    [ Info: 1a

"""
struct WebSocketLogger <: AbstractLogger
    stream::IO
    min_level::LogLevel
    meta_formatter
    show_limited::Bool
    right_justify::Int
    message_limits::Dict{Any,Int}
    shouldlog::Function
end
function WebSocketLogger(stream::IO=stderr, min_level=Info;
                       meta_formatter = default_metaformat, show_limited=true,
                       right_justify=0, shouldlog = shouldlog_default)
    WebSocketLogger(stream, min_level, meta_formatter,
                  show_limited, right_justify, Dict{Any,Int}(), shouldlog)
end

"Defines a default logging message format with timestamp"
function default_metaformat(level, _module, group, id, file, line)
    color = default_logcolor(level)
    if level == Wslog
        # Message generated with @wslog
        prefix = "Wslog $(Dates.Time(now())):"
        return color, prefix, ""
    elseif level == Info
        return color, string(level) * ':', ""
    else
        prefix = string(level) * ':'
        suffix = "$(Dates.Time(now())) @ "
        _module !== nothing && (suffix *= "$(_module)")
        if file !== nothing
            _module !== nothing && (suffix *= " ")
            suffix *= Base.contractuser(file)
            if line !== nothing
                suffix *= ":$(isa(line, UnitRange) ? "$(first(line))-$(last(line))" : line)"
            end
        end
        return color, prefix, suffix
    end
end

"Redirects to function that can be passed in when defining logger"
function shouldlog(logger::WebSocketLogger, level, _module, group, id)
    logger.shouldlog(logger::WebSocketLogger, level, _module, group, id)
end

"Early filtering of messages based on message id limits, silencing of HTTP.Servers if defined"
function shouldlog_default(logger::WebSocketLogger, level, _module, group, id)
    _module == WebSockets.HTTP.Servers && return false
    return get(logger.message_limits, id, 1) > 0
end

"Immutable lower bound for log level of accepted events."
min_enabled_level(logger::WebSocketLogger) = logger.min_level

"This is a copy of stdlib / Logging.ConsoleLogger's handle message function,
with the exceptions

    1) adding  :wslog=>true to the IOContext
    2) passing IOContext while converting message to strings
    3) evaluate Tuple message arguments without showing paranthesis and commas "
function handle_message(logger::WebSocketLogger, level, message, _module, group, id,
                        filepath, line; maxlog=nothing, kwargs...)
    if maxlog != nothing && maxlog isa Integer
        remaining = get!(logger.message_limits, id, maxlog)
        logger.message_limits[id] = remaining - 1
        remaining > 0 || return
    end

    # Generate a text representation of the message and all key value pairs,
    # split into lines.
    msgpretty = string_with_env_ws(logger.stream, message)
    msglines = [(indent=0,msg=l) for l in split(chomp(msgpretty), '\n')]
    dsize = displaysize(logger.stream)
    if !isempty(kwargs)
        valbuf = IOBuffer()
        rows_per_value = max(1, dsize[1]÷(length(kwargs)+1))
        valio = IOContext(IOContext(valbuf, logger.stream),
                          :displaysize => (rows_per_value,dsize[2]-5),
                          :limit => logger.show_limited,
                          :wslog=>true)
        for (key,val) in pairs(kwargs)
            showvalue(valio, val)
            vallines = split(String(take!(valbuf)), '\n')
            if length(vallines) == 1
                push!(msglines, (indent=2,msg=SubString("$key = $(vallines[1])")))
            else
                push!(msglines, (indent=2,msg=SubString("$key =")))
                append!(msglines, ((indent=3,msg=line) for line in vallines))
            end
        end
    end

    # Format lines as text with appropriate indentation and with a box
    # decoration on the left.
    color,prefix,suffix = logger.meta_formatter(level, _module, group, id, filepath, line)
    minsuffixpad = 2
    buf = IOBuffer()
    iob = IOContext(buf, logger.stream)
    nonpadwidth = 2 + (isempty(prefix) || length(msglines) > 1 ? 0 : length(prefix)+1) +
                  msglines[end].indent + termlength(msglines[end].msg) +
                  (isempty(suffix) ? 0 : length(suffix)+minsuffixpad)
    justify_width = min(logger.right_justify, dsize[2])
    if nonpadwidth > justify_width && !isempty(suffix)
        push!(msglines, (indent=0,msg=SubString("")))
        minsuffixpad = 0
        nonpadwidth = 2 + length(suffix)
    end
    for (i,(indent,msg)) in enumerate(msglines)
        boxstr = length(msglines) == 1 ? "[ " :
                 i == 1                ? "┌ " :
                 i < length(msglines)  ? "│ " :
                                         "└ "
        printstyled(iob, boxstr, bold=true, color=color)
        if i == 1 && !isempty(prefix)
            printstyled(iob, prefix, " ", bold=true, color=color)
        end
        print(iob, " "^indent, msg)
        if i == length(msglines) && !isempty(suffix)
            npad = max(0, justify_width - nonpadwidth) + minsuffixpad
            print(iob, " "^npad)
            printstyled(iob, suffix, color=:light_black)
        end
        println(iob)
    end

    write(logger.stream, take!(buf))
    nothing
end

"""
Context-aware text representation of the first argument to logging macros.
Made easily available for overloading on specific types.
"""
function string_with_env_ws(env, xs::Tuple)
    s = ""
    if isempty(xs)
        return s
    end
    for x in xs
        s *= string_with_env_ws(env, x)
    end
    return s
end
string_with_env_ws(env, x) = string_with_env(env, x)

macro  wslog(exs...) logmsg_code((@_sourceinfo)..., Wslog,  exs...) end


text:
function sel(n::Int64)
    return list -> list[n]
end

S1 = sel(1)
S2 = sel(2)
S3 = sel(3)

S1([1,2,3])

function larBernsteinBasis(u::Function, n::Int64)
    function map_fn(x::Array{Float64,1})
        bernstein(i, n, t::Float64) = binomial(n,i)*(1-t)^(n-i)*t^i
        return [bernstein(i, n, u(x)) for i in 0:n]
    end
    return map_fn
end

function larTensorProdSurface{F<:Function, T<:Real}(args::Array{F, 1},
                            controlpoints::Array{Array{Array{T, 1}, 1}, 1})
    ubasis , vbasis = args
    function map_fn(point::Array{Float64, 1})
        u, v = point
        U = ubasis([u])
        V = vbasis([v])
        target_dim = length(controlpoints[1][1])
        ret = zeros(Float64, target_dim)
        dim_u = length(U)
        dim_v = length(V)
        if (length(controlpoints)!= dim_u || length(controlpoints[1])!= dim_v)
            error("Invalid set of control points for those bases!")
        end
        for i in 1:dim_u
            for j in 1:dim_v
                #for M in 1:target_dim
                    for M in 1:target_dim
                        ret[M] += U[i]*V[j]*controlpoints[i][j][M]
                    end
                #end
            end
        end
        return ret
    end
    return map_fn
end

function larBilinearSurface(controlpoints)
   basis = larBernsteinBasis(sel(1), 1)
   return larTensorProdSurface([basis, basis], controlpoints)
end

function larBiquadraticSurface(controlpoints)
   basis1 = larBernsteinBasis(sel(1), 2)
   basis2 = larBernsteinBasis(sel(1), 2)
   return larTensorProdSurface([basis1,basis2], controlpoints)
end

function larBicubicSurface(controlpoints)
   basis1 = larBernsteinBasis(sel(1), 3)
   basis2 = larBernsteinBasis(sel(1), 3)
   return larTensorProdSurface([basis1, basis2], controlpoints)
end

function larBezier{T<:Real}(u::Function, controldata::Array{Array{T,1},1})
    n = length(controldata) - 1
    dim = length(controldata[1])
    function map_fn(point::Array{Float64, 1})
        t = u(point)
        out = zeros(Float64, dim)
        for i in 0:n
            weight = binomial(n,i)*(1-t)^(n-i)*t^i
            for k in 1:dim
                out[k] += weight*controldata[i+1][k]
            end
        end
        return out
     end
    return map_fn
end

function larBezierCurve(controlpoints)
   return larBezier(sel(1), controlpoints)
end

function larCoonsPatch{T<:Function}(args::Array{T, 1})
    su0_fn , su1_fn , s0v_fn , s1v_fn = args
    function map_fn(point)
        u, v = point
        su0 = su0_fn(point)
        su1 = su1_fn(point)
        s0v = s0v_fn(point)
        s1v = s1v_fn(point)
        ret = zeros(Float64, length(su0))
        for k in 1:length(ret)
            ret[k] = ((1-u)*s0v[k] + u*s1v[k]+(1-v)*su0[k] + v*su1[k] +
                     (1-u)*(1-v)*s0v[k] + (1-u)*v*s0v[k] + u*(1-v)*s1v[k] + u*v*s1v[k])
        end
        return ret
    end
    return map_fn
end

function bsplineBasis{T<:Real}(degree::Int64, knots::Array{T,1}, ncontrols::Int64)
    n = ncontrols - 1
    m = length(knots) - 1
    k = degree + 1
    tmin, tmax = knots[k-1], knots[n+1]
    if length(knots)!=(n+k+1)
        error("Invalid point/knots/degree for bspline!")
    end
    # de Boor coefficients
    function deBoor(i::Int64, k::Int64, t::Float64)
        # deBoor coefficient Ni1(t)
        if k == 1
            if(t >= knots[i] && t<knots[i+1]) || (t == tmax && t>=knots[i] && t<=knots[i+1])
                # t in [T_i, T_(i+1)) vel t=T_(n+1) in [T_i, T_(i+1)]
                return 1
            else
                return 0
            end
        end
        # deBoor coefficient Nik(t)
        ret = 0
        num1, div1 = t-knots[i], knots[i+k-1]-knots[i]
        if div1 != 0
            ret += (num1/div1) * deBoor(i,k-1,t)
        end
        num2, div2 = knots[i+k]-t, knots[i+k]-knots[i+1]
        if div2 != 0
            ret += (num2/div2) * deBoor(i+1,k-1,t)
        end
        return ret
    end
    # map function
    function map_fn{T<:Real}(point::Array{T, 1})
        t = point[1]
        return [deBoor(i,k,t) for i in 1:n+1]
    end
    return map_fn
end

function t_bspline{T<:Real, F<:Function}(u::Function, degree::Int64, knots::Array{T,1},
                                    points_fn::Array{F, 1})

    n = length(points_fn)-1
    m = length(knots)-1
    k = degree + 1
    tmin,tmax = knots[k-1], knots[n+1]

    # see http://www.na.iac.cnr.it/~bdv/cagd/spline/B-spline/bspline-curve.html
    if length(knots)!=(n+k+1)
        error("Invalid point/knots/degree for bspline!")
    end
    # de Boor coefficients
    function deBoor(i::Int64, k::Int64, t::Float64)
        # Ni1(t)
        if k == 1
            if(t>=knots[i] && t<knots[i+1]) || (t==tmax && t>=knots[i] && t<=knots[i+1])
                # i use strict inclusion for the max value
                return 1
            else
                return 0
            end
        end
        # Nik(t)
        ret = 0

        num1, div1 = t-knots[i], knots[i+k-1]-knots[i]
        if div1 != 0
            ret+=(num1/div1) * deBoor(i,k-1,t)
            # elif num1!=0: ret+=deBoor(i,k-1,t)
        end

        num2, div2 = knots[i+k]-t, knots[i+k]-knots[i+1]
        if div2!=0
            ret+=(num2/div2) * deBoor(i+1,k-1,t)
            # elif num2!=0: ret+=deBoor(i,k-1,t)
        end

        return ret
    end

    # map function
    function map_fn(point::Array{Float64, 1})
        t = u(point)
        # if control points are functions
        #points = point -> [f(point) for f in points_fn]    #POSSIBILE ERRORE QUI
        target_dim = length(points_fn[1](point))    #errore qui
        ret = zeros(Float64, target_dim)
        for i in 1:n+1
            coeff = deBoor(i,k,t)
            for m in 1:target_dim
                ret[m] += points_fn[i](point)[m] * coeff
            end
        end
        return ret
    end
    return map_fn
end

function t_nurbs{T<:Real, F<:Function}(u::Function, degree::Int64, knots::Array{T,1},
                                       points::Array{F, 1})
    bspline = t_bspline(u, degree, knots, points)
    function map_fn(point)          #eliminare
        ret = bspline(point)
        last = ret[end]  #last element in ret
        if last != 0
            ret = [value/last for value in ret]
        end
        ret = ret[1:end-1] #it slices off the last element (that is now 1)
        return ret
    end
    return map_fn
end


text:
import Base: getindex, setindex!, permutedims, permutedims!, size, copy, !, show, +, inv, similar

function check_species_validity(::Type{T}) where {T <: Any}
  throw(ArgumentError("The type $(T) is not an allowed species type"))
end

function check_species_validity(::Type{T}) where {T <: Union{Symbol,String}}
end

"""
    show(io::IO, N::AbstractEcologicalNetwork)
"""
function show(io::IO, N::AbstractEcologicalNetwork)
  p = typeof(N) <: AbstractBipartiteNetwork ? "bipartite" : "unipartite"
  t = ""
  t = typeof(N) <: ProbabilisticNetwork ? "probabilistic" : t
  t = typeof(N) <: QuantitativeNetwork ? "quantitative" : t
  print(io, "$(richness(N; dims=1))×$(richness(N; dims=2)) $(p) $(t) ecological network $(eltype(N)) (L: $(links(N)))")
end

"""
    richness(N::AbstractEcologicalNetwork, i::Int64)

Returns the number of species on either side of the network. The value of `i`
can be `1` (top-level species) or `2` (bottom-level species), as in the
`species` function.
"""
function richness(N::AbstractEcologicalNetwork; dims::Union{Nothing,Integer}=nothing)
  return length(species(N; dims=dims))
end

"""
    species(N::AbstractBipartiteNetwork)

Returns an array of all species in a bipartite network. The order of the species
corresponds to the order of rows (top level) and columns (bottom level) of the
adjacency matrix, in this order.
"""
function species(N::AbstractBipartiteNetwork; dims::Union{Nothing,Integer}=nothing)
  dims === nothing && return vcat(N.T, N.B)
  dims == 1 && return N.T
  dims == 2 && return N.B
  throw(ArgumentError("dims can only be 1 (top species) or 2 (bottom species), or `nothing` (all species), you used $(dims)"))
end

"""
    species(N::AbstractUnipartiteNetwork; dims::Int64=1)

Returns an array of species on either side of a unipartite network. In a
unipartite network, species are the same on either size. This function is
nevertheless useful when you want to write code that takes either side of the
network in a general way.
"""
function species(N::AbstractUnipartiteNetwork; dims::Union{Nothing,Integer}=nothing)
  return N.S
end

function species_objects(N::AbstractUnipartiteNetwork)
  return (N.S,)
end

function species_objects(N::AbstractBipartiteNetwork)
  return (N.T, N.B)
end

"""
    has_interaction{(N::AbstractEcologicalNetwork, i::NT, j::NT)

This function returns `true` if the interaction between `i` and `j` is not 0.
This refers to species by their names/values, and is the recommended way to test
for the presence of an interaction.

Use `N[i,j]` if you need to get the value of the interaction.
"""
function has_interaction(N::AbstractEcologicalNetwork, i::NT, j::NT) where {NT}
  check_species_validity(NT)
  @assert i ∈ species(N; dims=1)
  @assert j ∈ species(N; dims=2)
  i_pos = something(findfirst(isequal(i), species(N; dims=1)),0)
  j_pos = something(findfirst(isequal(j), species(N; dims=2)),0)
  return has_interaction(N, i_pos, j_pos)
end

"""
    has_interaction(N::AbstractEcologicalNetwork, i::Int64, j::Int64)

This function returns `true` if the interaction between `i` and `j` is not 0.
This refers to species by their position instead of their name, and is not
recommended as the main solution. This is used internally by a few functions,
but is exported because it may be of general use.
"""
function has_interaction(N::AbstractEcologicalNetwork, i::Int64, j::Int64)
  # We need to make sure that the interaction is somewhere within the network
  @assert i <= size(N.A, 1)
  @assert j <= size(N.A, 2)
  # This should be reasonably general...
  return N[i,j] != zero(typeof(N[i,j]))
end

"""
    nodiagonal!(N::AbstractUnipartiteNetwork)

Modifies the network so that its diagonal is set to the appropriate zero.
"""
function nodiagonal!(N::AbstractUnipartiteNetwork)
  for i in 1:richness(N)
    N.A[i,i] = zero(eltype(N.A))
  end
end

"""
    nodiagonal(N::AbstractUnipartiteNetwork)

Returns a copy of the network with its diagonal set to zero.
"""
function nodiagonal(N::AbstractUnipartiteNetwork)
  Y = copy(N)
  nodiagonal!(Y)
  return Y
end

"""
    nodiagonal(N::AbstractBipartiteNetwork)

Returns a *copy* of the network (because the diagonal of a bipartite network is
never a meaningful notion). This function is clearly useless, but allows to
write general code for all networks types when a step requires removing the
diagonal.
"""
function nodiagonal(N::AbstractBipartiteNetwork)
  return copy(N)
end

"""
    nodiagonal!(N::AbstractBipartiteNetwork)

Does nothing.
"""
function nodiagonal!(N::AbstractBipartiteNetwork)
end

"""
    size(N::AbstractEcologicalNetwork)

Return the size of the adjacency matrix of an `AbstractEcologicalNetwork` object.
"""
function size(N::AbstractEcologicalNetwork)
  size(N.A)
end

"""
Creates a copy of a network -- this returns an object with the same type, and
the same content.
"""
function copy(N::AbstractEcologicalNetwork)
  a = copy(N.A)
  sp = copy.(species_objects(N))
  return typeof(N)(a, sp...)
end

"""
    getindex(N::AbstractEcologicalNetwork, i::T, j::T)

Get the value of an interaction based on the position of the species.
"""
function getindex(N::AbstractEcologicalNetwork, i::T, j::T) where {T <: Int}
  return N.A[i,j]
end

"""
    getindex(N::AbstractEcologicalNetwork, i::T)

Get the value of an interaction based on the position of the species.
"""
function getindex(N::AbstractEcologicalNetwork, i::T) where {T <: Int}
  return N.A[i]
end

"""
    getindex(N::AbstractEcologicalNetwork, ::Colon, j::T)

Get the value of an interaction based on the position of the species.
"""
function getindex(N::AbstractEcologicalNetwork, ::Colon, j::T) where {T <: Int}
  return N.A[:,j]
end

"""
    getindex(N::AbstractEcologicalNetwork, i::T, ::Colon)

Get the value of an interaction based on the position of the species.
"""
function getindex(N::AbstractEcologicalNetwork, i::T, ::Colon) where {T <: Int}
  return N.A[i,:]
end

"""
    getindex{T}(N::AbstractEcologicalNetwork, s1::T, s2::T)

Get the value of an interaction based on the *name* of the species. This is the
recommended way to look for things in a network.
"""
function getindex(N::AbstractEcologicalNetwork, s1::T, s2::T) where {T}
  @assert T == last(eltype(N))
  @assert s1 ∈ species(N; dims=1)
  @assert s2 ∈ species(N; dims=2)
  s1_pos = something(findfirst(isequal(s1), species(N; dims=1)), 0)
  s2_pos = something(findfirst(isequal(s2), species(N; dims=2)), 0)
  return N[s1_pos, s2_pos]
end

"""
    getindex{T}(N::AbstractBipartiteNetwork, ::Colon, sp::T)

Gets the predecessors (*i.e.* species that interacts with / consume) of a focal
species. This returns the list of species as a `Set` object, in which ordering
is unimportant.
"""
function getindex(N::AbstractEcologicalNetwork, ::Colon, sp::T) where {T}
  @assert T == last(eltype(N))
  @assert sp ∈ species(N; dims=2)
  return Set(filter(x -> has_interaction(N, x, sp), species(N; dims=1)))
end

"""
    getindex{T}(N::AbstractEcologicalNetwork, sp::T, ::Colon)

Gets the successors (*i.e.* species that are interacted with / consumed) of a
focal species. This returns the list of species as a `Set` object, in which
ordering is unimportant.
"""
function getindex(N::AbstractEcologicalNetwork, sp::T, ::Colon) where {T}
  @assert T == last(eltype(N))
  @assert sp ∈ species(N; dims=1)
  return Set(filter(x -> has_interaction(N, sp, x), species(N; dims=2)))
end

"""
    getindex{T}(N::AbstractUnipartiteNetwork, sp::Array{T})

Induce a unipartite network based on a list of species, all of which must be in
the original network. This function takes a single argument (as opposed to two
arrays, or an array and a colon) to ensure that the returned network is
unipartite.

The network which is returned by this function may not have the species in the
order specified by the user for performance reasons.
"""
function getindex(N::AbstractUnipartiteNetwork, sp::Vector{T}) where {T}
  @assert T == last(eltype(N))
  @assert all(map(x -> x ∈ species(N), sp))
  sp_pos = findall((in)(sp), species(N))
  n_sp = species(N)[sp_pos]
  n_int = N.A[sp_pos, sp_pos]
  return typeof(N)(n_int, n_sp)
end

"""
    getindex{T}(N::AbstractBipartiteNetwork, ::Colon, sp::Array{T})

TODO
"""
function getindex(N::AbstractBipartiteNetwork, ::Colon, sp::Vector{T}) where {T}
  @assert T == last(eltype(N))
  @assert all(map(x -> x ∈ species(N; dims=2), sp))
  sp_pos = findall((in)(sp), species(N; dims=2))
  n_t = N.T
  n_b = N.B[sp_pos]
  n_int = N.A[:, sp_pos]
  return typeof(N)(n_int, n_t, n_b)
end

"""
    getindex{T}(N::AbstractBipartiteNetwork, sp::Array{T}, ::Colon)

TODO
"""
function getindex(N::AbstractBipartiteNetwork, sp::Vector{T}, ::Colon) where {T}
  @assert T == last(eltype(N))
  @assert all(map(x -> x ∈ species(N; dims=1), sp))
  sp_pos = findall((in)(sp), species(N; dims=1))
  n_t = N.T[sp_pos]
  n_b = N.B
  n_int = N.A[sp_pos,:]
  return typeof(N)(n_int, n_t, n_b)
end

"""
    getindex{T}(N::AbstractBipartiteNetwork, sp1::Array{T}, sp2::Array{T})

TODO
"""
function getindex(N::AbstractBipartiteNetwork, sp1::Vector{T}, sp2::Vector{T}) where {T <: Integer}
   @assert all(sp1 .<= richness(N; dims=1))
   @assert all(sp2 .<= richness(N; dims=2))
   n_t = N.T[sp1]
   n_b = N.B[sp2]
   n_int = N.A[sp1,sp2]
   return typeof(N)(n_int, n_t, n_b)
end

"""
    getindex{T}(N::AbstractBipartiteNetwork, sp1::Array{T}, sp2::Array{T})

TODO
"""
function getindex(N::AbstractBipartiteNetwork, sp1::Vector{T}, sp2::Vector{T}) where {T}
  @assert T == last(eltype(N))
  @assert all(map(x -> x ∈ species(N; dims=1), sp1))
  @assert all(map(x -> x ∈ species(N; dims=2), sp2))
  sp1_pos = findall((in)(sp1), species(N; dims=1))
  sp2_pos = findall((in)(sp2), species(N; dims=2))
  return N[sp1_pos, sp2_pos]
end

"""
    setindex!(N::T, A::Any, i::E, j::E) where {T <: AbstractEcologicalNetwork, E}

Changes the value of the interaction at the specificied position, where `i` and
`j` are species *names*. Note that this operation **changes the network**.
"""
function setindex!(N::T, A::Any, i::E, j::E) where {T <: AbstractEcologicalNetwork, E}
  @assert E == last(eltype(N))
  @assert i ∈ species(N; dims=1)
  @assert j ∈ species(N; dims=2)
  i_pos = something(findfirst(isequal(i), species(N; dims=1)),0)
  j_pos = something(findfirst(isequal(j), species(N; dims=2)),0)
  N.A[i_pos, j_pos] = A
end

"""
    setindex!(N::T, A::K, i::E, j::E) where {T <: AbstractEcologicalNetwork, K <: first(eltype(N)), E <: Int}

Changes the value of the interaction at the specificied position, where `i` and
`j` are species *positions*. Note that this operation **changes the network**.
"""
function setindex!(N::T, A::Any, i::E, j::E) where {T <: AbstractEcologicalNetwork, E <: Int}
  @assert i ≤ richness(N; dims=1)
  @assert j ≤ richness(N; dims=2)
  N.A[i, j] = A
end

"""
    Base.:!{T<:DeterministicNetwork}(N::T)

Returns the inverse of a binary network -- interactions that were `false` become
`true`, and conversely.
"""
function Base.:!(N::T) where {T <: BinaryNetwork}
  return typeof(N)(.!N.A, species_objects(N)...)
end

"""
    permutedims(N::AbstractBipartiteNetwork)

Tranposes the network and returns a copy
"""
function permutedims(N::AbstractEcologicalNetwork)
  NA = permutedims(N.A)
  new_sp = typeof(N) <: AbstractBipartiteNetwork ? (N.B, N.T) : (N.S,)
  return typeof(N)(NA, new_sp...)
end

"""
    inv(N::AbstractUnipartiteNetwork)

Inverses interactions in a network. Note that this is *not* the same as
changing the side of species (`permutedims`); this function is mostly useful
when the fluxes of biomass are represented as going from the resource to
the consumer, which is the opposite direction of the interaction.
"""
function inv(N::AbstractUnipartiteNetwork)
   Y = copy(N)
   inv!(Y)
   return Y
end

"""
    inv!(N::AbstractUnipartiteNetwork)

Inverses interactions in a network. Note that this is *not* the same as
changing the side of species (`permutedims`); this function is mostly useful
when the fluxes of biomass are represented as going from the resource to the
consumer, which is the opposite direction of the interaction. This functions
modifies the network.
"""
function inv!(N::AbstractUnipartiteNetwork)
   for i in 1:(richness(N)-1)
      for j in (i+1):richness(N)
         N.A[i,j] = N.A[j,i]
      end
   end
end

"""
    interactions(N::AbstractEcologicalNetwork)

Returns the interactions in the ecological network. Interactions are returned as
an array of named tuples. *A minima*, these have fields `from` and `to`. For
networks that are probabilistic, there is a `probability` field. For networks
that are quantitative, there is a `strength` field. This functions allows to
iterate over interactions in a network in a convenient way.
"""
function interactions(N::AbstractEcologicalNetwork)
  edges_accumulator = NamedTuple[]
  fields = [:from, :to]
  if typeof(N) <: ProbabilisticNetwork
    push!(fields, :probability)
  end
  if typeof(N) <: QuantitativeNetwork
    push!(fields, :strength)
  end
  sp1 = species(N; dims=1)
  sp2 = species(N; dims=2)
  for i in eachindex(sp1)
    s1 = sp1[i]
    for j in eachindex(sp2)
      if has_interaction(N, i, j)
        values = Any[s1, sp2[j]]
        if typeof(N) <: ProbabilisticNetwork
          push!(values, N[i,j])
        end
        if typeof(N) <: QuantitativeNetwork
          push!(values, N[i,j])
        end
        int_nt = NamedTuple{tuple(fields...)}(tuple(values...))
        push!(edges_accumulator, int_nt)
      end
    end
  end
  return unique(edges_accumulator)
end

"""
    +(n1::T, n2::T) where {T <: BipartiteQuantitativeNetwork}

Adds two quantitative bipartite networks. TODO
"""
function +(n1::T, n2::T) where {T <: BipartiteQuantitativeNetwork}
  st = union(species(n1; dims=1), species(n2; dims=1))
  sb = union(species(n1; dims=2), species(n2; dims=2))
  A = zeros(first(eltype(n1)), (length(st), length(sb)))
  N = T(A, st, sb)
  for i1 in n1
    N[i1.from,i1.to] = N[i1.from,i1.to] + i1.strength
  end
  for i2 in n2
    N[i2.from,i2.to] = N[i2.from,i2.to] + i2.strength
  end
  return N
end

"""
    similar(N::T) where {T <: AbstractEcologicalNetwork}

Returns a network with the same species, but an empty interaction matrix. This
is useful if you want to generate a "blank slate" for some analyses.
"""
function similar(N::T) where {T <: AbstractEcologicalNetwork}
   Y = zeros(first(eltype(N)), size(N))
   return T(Y, species_objects(N)...)
end


text:
using Printf
import LinearAlgebra: ×

export AbstractDomain
export AbstractBoundary

export DomainBoundary
export IntervalDomain, ProductDomain, Circle, S¹, Interval
export ×, dim, info, check_full_periodicity

abstract type AbstractDomain end
abstract type AbstractBoundary end

struct DomainBoundary <: AbstractBoundary
    closure
end
struct IntervalDomain{AT, BT, PT} <: AbstractDomain
    a::AT
    b::BT
    periodic::PT
end

function IntervalDomain(a, b; periodic=false)
    @assert a < b
    return IntervalDomain(a, b, periodic)
end
function Circle(a, b)
    @assert a < b
    return IntervalDomain(a, b, periodic = true)
end
S¹ = Circle
function Interval(a, b)
    @assert a < b
    return IntervalDomain(a, b)
end

function Base.show(io::IO, Ω::IntervalDomain)
    a = Ω.a
    b = Ω.b
    printstyled(io, "[", color = 226)
    printstyled("$a, $b", color = 7)
    Ω.periodic ? printstyled(io, ")", color = 226) : printstyled(io, "]", color = 226)
 end

∂(Ω::IntervalDomain) = Ω.periodic ? DomainBoundary(nothing) : DomainBoundary((Ω.a, Ω.b))

struct ProductDomain{DT} <: AbstractDomain
    domains::DT
end

function Base.show(io::IO, Ω::ProductDomain)
    for (i,domain) in enumerate(Ω.domains)
        print(domain)
        if i != length(Ω.domains)
            printstyled(io, "×", color = 118)
        end
    end
 end

 function dim(Ω::IntervalDomain)
    return 1
end

function dim(Ω::ProductDomain)
    return length(Ω.domains)
end

×(arg1::AbstractDomain, arg2::AbstractDomain) = ProductDomain((arg1, arg2))
×(args::ProductDomain, arg2::AbstractDomain) = ProductDomain((args.domains..., arg2))
×(arg1::AbstractDomain, args::ProductDomain) = ProductDomain((arg1, args.domains...))
×(arg1::ProductDomain, args::ProductDomain) = ProductDomain((arg1.domains..., args.domains...))
×(args::AbstractDomain) = ProductDomain(args...)

function info(Ω::ProductDomain)
    println("This is a ", dim(Ω),"-dimensional tensor product domain.")
    print("The domain is ")
    println(Ω, ".")
    for (i,domain) in enumerate(Ω.domains)
        domain_string = domain.periodic ? "periodic" : "wall-bounded"
        length = @sprintf("%.2f ", domain.b-domain.a)
        println("The dimension $i domain is ", domain_string, " with length ≈ ", length)

    end
    return nothing
end

function check_full_periodicity(Ω::ProductDomain)
    b = [Ω.domains[i].periodic for i in eachindex(Ω.domains)]
    return prod(b)
end

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
prompt_1 = ""*
"This tutorial will explain to solve linear regression using the ordinary least squares (OLS) in Julia."*
"First, we need dataset with x and y variables to run regression with one variable."*

request_1 = "Create two variables, x and y. To make this simple, just copy paste this code: x = [10,12,11,6,7,23,14,15]"
evaluator_1 = ExpressionEvaluator(:x, [10,12,11,6,7,23,14,15])

request_2 = "Now copy paste this: y = [234,521,311,201,102,71,111,900]"
evaluator_2 = ExpressionEvaluator(:y, [234,521,311,201,102,71,111,900])

prompt_2 = ""*
"Please add both x and y in one data frame. Install and load the DataFrames.jl package, as we did for GLM.jl and then exectue the following code to add them in one data frame named `data`"*
"`data = DataFrame(x=x, y=y)`"*

request_3 = ""
evaluator_3 =

prompt_3 = ""*
"Now we need GLM.jl package for solving regression model using OLS"*
"Execute:"*
"\n"*
"using Pkg; Pkg.add(GLM); using GLM"*
"GLM should be in between uottes"*

request_4 = "For solving this regression, we can use `lm()` from GLM.jl package. Run: ols = lm(@formula(y~x), data)"
evaluator_4 = ExpressionEvaluator(:ols, ols)

lesson_1 = Lesson( prompt_1, request_1, evaluator_1, request_2, evaluator_2, prompt_2, request_3, evaluator_3, prompt_3, request_4, evaluator_4)

lesson_plan = [ lesson_1]


text:
module OscarDevTools

import GitHub
import JSON
import TOML
import Pkg
import LibGit2

export oscar_develop, oscar_update, oscar_branch, oscar_add_remotes

include("helpers.jl")
include("githubci.jl")
include("devutils.jl")

end # module


text:
# Eric S. Tellez <eric.tellez@infotec.mx>

#Pkg.add("GZip")
#Pkg.add("Glob")
#Pkg.add("JSON")

using ArgParse
import JSON
import Glob

include("textmodel.jl")
include("io.jl")

function groupby(fun, input, textkey, userkey, popkey, minsize)
    H = Dict{String,Dict}()
    lineno = 0

    dopopkey = popkey != ""
    itertweets(input) do tweet
        lineno += 1
        key = getkeypath(tweet, userkey)
        text = getkeypath(tweet, textkey)

        if dopopkey
            arr = split(key, popkey[1])
            pop!(arr)
            key = join(arr, popkey[1])
        end

        if haskey(H, key)
            push!(H[key][textkey], text)
        else
            tweet[textkey] = [text]
            tweet[userkey] = key
            # info(JSON.json(tweet), " --- ", length(H))
            H[key] = tweet
        end

        if lineno % 1000 == 0
            info("Advance $uniques from $lineno")
        end
    end

    for (k, v) in H
        if length(v[textkey]) >= minsize
            fun(v)
        end
    end
end

function main()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--by"
        help = "specifies the key to group"
        arg_type = String
        default = "key"
        "--text"
        help = "specifies the key pointing to the text"
        arg_type = String
        default = "text"
        "--pop"
        help = "specifies the character to split the key; also it applies a pop to create the key"
        arg_type = String
        default = ""
        "--minsize"
        help = "specifies the minimum size of the group to be part of the output"
        arg_type = Int
        default = 1

    end

    args = parse_args(ARGS, s)
    groupby(STDIN, args["text"], args["by"], args["pop"], args["minsize"]) do tweet
        write(STDOUT, JSON.json(tweet), '\n')
    end
end

main()

max_stars_repo_name: INGEOTEC/PatternAnalysisTools

max_issues_repo_name: INGEOTEC/PatternAnalysisTools

max_forks_repo_name: INGEOTEC/PatternAnalysisTools


text:
#Ax_mul_Bx

export Ax_mul_Bx

"""
`Ax_mul_Bx(A::AbstractOperator,B::AbstractOperator)`

Create an operator `P` such that:

`P*x == (Ax)*(Bx)`

# Example

```julia
julia> A,B = randn(4,4),randn(4,4);

julia> P = Ax_mul_Bx(MatrixOp(A,4),MatrixOp(B,4))
▒*▒  ℝ^4 -> ℝ^(4, 4)

julia> X = randn(4,4);

julia> P*X == (A*X)*(B*X)
true

```
"""
struct Ax_mul_Bx{
                 L1 <: AbstractOperator,
                 L2 <: AbstractOperator,
                 C <: AbstractArray,
                 D <: AbstractArray,
                } <: NonLinearOperator
  A::L1
  B::L2
  bufA::C
  bufB::C
  bufC::C
  bufD::D
  function Ax_mul_Bx(A::L1, B::L2, bufA::C, bufB::C, bufC::C, bufD::D) where {L1,L2,C,D}
    if ndims(A,1) != 2 || size(A,2) != size(B,2) || size(A,1)[2] != size(B,1)[1]
      throw(DimensionMismatch("Cannot compose operators"))
    end
    new{L1,L2,C,D}(A,B,bufA,bufB,bufC,bufD)
  end
end

struct Ax_mul_BxJac{
                    L1 <: AbstractOperator,
                    L2 <: AbstractOperator,
                    C <: AbstractArray,
                    D <: AbstractArray,
                   } <: LinearOperator
  A::L1
  B::L2
  bufA::C
  bufB::C
  bufC::C
  bufD::D
end

# Constructors
function Ax_mul_Bx(A::AbstractOperator,B::AbstractOperator)
  s,t = size(A,1), codomainType(A)
  bufA = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...)
  s,t = size(B,1), codomainType(B)
  bufB = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...)
  bufC = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...)
  s,t = size(A,2), domainType(A)
  bufD = eltype(s) <: Int ? zeros(t,s) : ArrayPartition(zeros.(t,s)...)
  Ax_mul_Bx(A,B,bufA,bufB,bufC,bufD)
end

# Jacobian
function Jacobian(P::Ax_mul_Bx{L1,L2,C,D}, x::AbstractArray) where {L1,L2,C,D}
  JA, JB = Jacobian(P.A, x), Jacobian(P.B, x)
  Ax_mul_BxJac{typeof(JA),typeof(JB),C,D}(JA,JB,P.bufA,P.bufB,P.bufC,P.bufD)
end

# Mappings
function mul!(y, P::Ax_mul_Bx{L1,L2,C,D}, b) where {L1,L2,C,D}
  mul!(P.bufA,P.A,b)
  mul!(P.bufB,P.B,b)
  mul!(y,P.bufA, P.bufB)
end

function mul!(y, J::AdjointOperator{Ax_mul_BxJac{L1,L2,C,D}}, b) where {L1,L2,C,D}
  #y .= J.A.B' * ( J.A.bufA'*b ) + J.A.A' * ( b*J.A.bufB' )
  mul!(J.A.bufC, J.A.bufA', b)
  mul!(y, J.A.B', J.A.bufC)
  mul!(J.A.bufA, b, J.A.bufB')
  mul!(J.A.bufD, J.A.A', J.A.bufA)
  y .+= J.A.bufD
  return y
end

size(P::Union{Ax_mul_Bx,Ax_mul_BxJac}) = ((size(P.A,1)[1],size(P.B,1)[2]),size(P.A,2))

fun_name(L::Union{Ax_mul_Bx,Ax_mul_BxJac}) = fun_name(L.A)*"*"*fun_name(L.B)

domainType(L::Union{Ax_mul_Bx,Ax_mul_BxJac})   = domainType(L.A)
codomainType(L::Union{Ax_mul_Bx,Ax_mul_BxJac}) = codomainType(L.A)

# utils
function permute(P::Ax_mul_Bx{L1,L2,C,D},
                 p::AbstractVector{Int}) where {L1,L2,C,D <:ArrayPartition}
  Ax_mul_Bx(permute(P.A,p),permute(P.B,p),P.bufA,P.bufB,P.bufC,ArrayPartition(P.bufD.x[p]) )
end

remove_displacement(P::Ax_mul_Bx) =
Ax_mul_Bx(remove_displacement(P.A), remove_displacement(P.B), P.bufA, P.bufB, P.bufC, P.bufD)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "NOMAD"
version = v"1.0.0"

# Collection of sources required to build NOMAD
sources = [
    "https://www.gerad.ca/nomad/Downloads/unix_linux/NOMAD.zip" =>
    "2134e64d5c8728054a797d0067eb49083d66ca64ed8ad7e75c34700b528aeff3",

]

# Bash recipe for building across all platforms
script = raw"""
export NOMAD_HOME=${WORKSPACE}/srcdir/nomad.3.9.1
export PATH=${NOMAD_HOME}/bin:$PATH
cd $NOMAD_HOME
./configure
make
rm -rf doc
rm -rf bin
rm -rf examples
rm -rf lib
rm -rf utils
rm -rf tools
rm -rf ext/sgtelib/bin
rm -rf ext/sgtelib/example
rm -rf ext/sgtelib/matlab_server
rm -rf ext/sgtelib/user_guide
cd builds/release/lib
rm libsgtelib.so
ln -s ../../../ext/sgtelib/lib/libsgtelib.so libsgtelib.so
cp -rf ${NOMAD_HOME} ${prefix}
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc)
]

# The products that we will ensure are always built
products(prefix) = Array{Product,1}()

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)