text:
# MathOptInterface constraints

##################################################
## Support constraints
MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.ZeroOne}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{MOI.Integer}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.SingleVariable}, ::Type{<:SS}) = true

MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.LessThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.EqualTo{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{MOI.Interval{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarQuadraticFunction{Float64}}, ::Type{MOI.LessThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarQuadraticFunction{Float64}}, ::Type{MOI.GreaterThan{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{MOI.ScalarQuadraticFunction{Float64}}, ::Type{MOI.EqualTo{Float64}}) = true
MOI.supports_constraint(::Optimizer, ::Type{<:SF}, ::Type{<:SS}) = true
MOI.supports_constraint(::Optimizer, ::Type{VAF}, ::Type{<:VLS}) = true
MOI.supports_constraint(::Optimizer, ::Type{VOV}, ::Type{<:VLS}) = true
MOI.supports_constraint(::Optimizer, ::Type{VAF}, ::Type{MOI.SecondOrderCone}) = true
MOI.supports_constraint(::Optimizer, ::Type{VOV}, ::Type{MOI.SecondOrderCone}) = true
MOI.supports_constraint(::Optimizer, ::Type{VOV}, ::Type{MOI.Complements}) = true
MOI.supports_constraint(::Optimizer, ::Type{VAF}, ::Type{MOI.Complements}) = true

##################################################
## Getters
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints) =
    KN_get_number_cons(model.inner)
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{MOI.SingleVariable, MOI.ZeroOne}) =
    model.number_zeroone_constraints
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{MOI.SingleVariable, MOI.Integer}) =
    model.number_integer_constraints
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{MOI.SingleVariable, S}) where S <: LS =
    sum(typeof.(collect(keys(model.constraint_mapping))) .== MOI.ConstraintIndex{MOI.SingleVariable, S})
# TODO: a bit hacky, but that should work for MOI Test.
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{VAF, MOI.Nonnegatives}) =
    sum(typeof.(collect(keys(model.constraint_mapping))) .== MOI.ConstraintIndex{VAF, MOI.Nonnegatives})
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{VAF, MOI.Nonpositives}) =
    sum(typeof.(collect(keys(model.constraint_mapping))) .== MOI.ConstraintIndex{VAF, MOI.Nonpositives})
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{VAF, MOI.Zeros}) =
    sum(typeof.(collect(keys(model.constraint_mapping))) .== MOI.ConstraintIndex{VAF, MOI.Zeros})
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{VAF, MOI.SecondOrderCone}) =
    sum(typeof.(collect(keys(model.constraint_mapping))) .== MOI.ConstraintIndex{VAF, MOI.SecondOrderCone})
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{VOV, T}) where T <: VLS =
    sum(typeof.(collect(keys(model.constraint_mapping))) .== MOI.ConstraintIndex{VOV, T})
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64}, S}) where S <: LS  =
    sum(typeof.(collect(keys(model.constraint_mapping))) .== MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, S})
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{MOI.ScalarQuadraticFunction{Float64}, S}) where S <: LS  =
    sum(typeof.(collect(keys(model.constraint_mapping))) .== MOI.ConstraintIndex{MOI.ScalarQuadraticFunction{Float64}, S})

##################################################
# Generic constraint definition
#--------------------------------------------------
# Bound constraint on variables.
function MOI.add_constraint(model::Optimizer,
                            v::MOI.SingleVariable, lt::MOI.LessThan{Float64})
    (model.number_solved >= 1) && throw(AddConstraintError())
    vi = v.variable
    check_inbounds(model, vi)
    if isnan(lt.upper)
        error("Invalid upper bound value $(lt.upper).")
    end
    if has_upper_bound(model, vi)
        error("Upper bound on variable $vi already exists.")
    end
    if is_fixed(model, vi)
        error("Variable $vi is fixed. Cannot also set upper bound.")
    end
    ub = check_value(lt.upper)
    model.variable_info[vi.value].has_upper_bound = true
    # By construction, MOI's indexing is the same as KNITRO's indexing.
    KN_set_var_upbnds(model.inner, vi.value - 1, ub)
    ci = MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}}(vi.value)
    model.constraint_mapping[ci] = convert(Cint, vi.value)
    return ci
end

function MOI.add_constraint(model::Optimizer,
                            v::MOI.SingleVariable, gt::MOI.GreaterThan{Float64})
    (model.number_solved >= 1) && throw(AddConstraintError())
    vi = v.variable
    check_inbounds(model, vi)
    if isnan(gt.lower)
        error("Invalid lower bound value $(gt.lower).")
    end
    if has_lower_bound(model, vi)
        error("Lower bound on variable $vi already exists.")
    end
    if is_fixed(model, vi)
        error("Variable $vi is fixed. Cannot also set lower bound.")
    end
    lb = check_value(gt.lower)
    model.variable_info[vi.value].has_lower_bound = true
    # We assume that MOI's indexing is the same as KNITRO's indexing.
    KN_set_var_lobnds(model.inner, vi.value - 1, lb)
    ci = MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}}(vi.value)
    model.constraint_mapping[ci] = convert(Cint, vi.value)
    return ci
end

function MOI.add_constraint(model::Optimizer,
                            v::MOI.SingleVariable, set::MOI.Interval{Float64})
    (model.number_solved >= 1) && throw(AddConstraintError())
    vi = v.variable
    check_inbounds(model, vi)
    if isnan(set.lower) || isnan(set.upper)
        error("Invalid lower bound value $(set.lower).")
    end
    if has_lower_bound(model, vi) || has_upper_bound(model, vi)
        error("Bounds on variable $vi already exists.")
    end
    if is_fixed(model, vi)
        error("Variable $vi is fixed. Cannot also set lower bound.")
    end
    lb = check_value(set.lower)
    ub = check_value(set.upper)
    model.variable_info[vi.value].has_lower_bound = true
    model.variable_info[vi.value].has_upper_bound = true
    # We assume that MOI's indexing is the same as KNITRO's indexing.
    KN_set_var_lobnds(model.inner, vi.value - 1, lb)
    KN_set_var_upbnds(model.inner, vi.value - 1, ub)
    ci = MOI.ConstraintIndex{MOI.SingleVariable, MOI.Interval{Float64}}(vi.value)
    model.constraint_mapping[ci] = convert(Cint, vi.value)
    return ci
end

function MOI.add_constraint(model::Optimizer,
                            v::MOI.SingleVariable, eq::MOI.EqualTo{Float64})
    (model.number_solved >= 1) && throw(AddConstraintError())
    vi = v.variable
    check_inbounds(model, vi)
    if isnan(eq.value)
        error("Invalid fixed value $(eq.value).")
    end
    if has_lower_bound(model, vi)
        error("Variable $vi has a lower bound. Cannot be fixed.")
    end
    if has_upper_bound(model, vi)
        error("Variable $vi has an upper bound. Cannot be fixed.")
    end
    if is_fixed(model, vi)
        error("Variable $vi is already fixed.")
    end
    eqv = check_value(eq.value)
    model.variable_info[vi.value].is_fixed = true
    # We assume that MOI's indexing is the same as KNITRO's indexing.
    KN_set_var_fxbnds(model.inner, vi.value - 1, eqv)
    ci = MOI.ConstraintIndex{MOI.SingleVariable, MOI.EqualTo{Float64}}(vi.value)
    model.constraint_mapping[ci] = convert(Cint, vi.value)
    return ci
end

function MOI.add_constraint(model::Optimizer,
                            func::MOI.ScalarAffineFunction{Float64}, set::LS)
    (model.number_solved >= 1) && throw(AddConstraintError())
    check_inbounds(model, func)
    # Add a single constraint in KNITRO.
    num_cons = KN_add_con(model.inner)
    # Add bound to constraint.
    if isa(set, MOI.LessThan{Float64})
        val = check_value(set.upper)
        KN_set_con_upbnd(model.inner, num_cons, val - func.constant)
    elseif isa(set, MOI.GreaterThan{Float64})
        val = check_value(set.lower)
        KN_set_con_lobnd(model.inner, num_cons, val - func.constant)
    elseif isa(set, MOI.EqualTo{Float64})
        val = check_value(set.value)
        KN_set_con_eqbnd(model.inner, num_cons, val - func.constant)
    end
    # Parse structure of constraint.
    indexvars, coefs = canonical_linear_reduction(func)
    KN_add_con_linear_struct(model.inner, num_cons, indexvars, coefs)
    # Add constraint to index.
    ci = MOI.ConstraintIndex{typeof(func), typeof(set)}(num_cons)
    model.constraint_mapping[ci] = num_cons
    return ci
end

function MOI.add_constraint(model::Optimizer,
                            func::MOI.ScalarAffineFunction{Float64}, set::MOI.Interval{Float64})
    (model.number_solved >= 1) && throw(AddConstraintError())
    check_inbounds(model, func)
    # Add a single constraint in KNITRO.
    num_cons = KN_add_con(model.inner)
    # Add upper bound.
    lb = check_value(set.lower)
    ub = check_value(set.upper)
    KN_set_con_lobnd(model.inner, num_cons, lb - func.constant)
    KN_set_con_upbnd(model.inner, num_cons, ub - func.constant)
    # Parse structure of constraint.
    indexvars, coefs = canonical_linear_reduction(func)
    KN_add_con_linear_struct(model.inner, num_cons, indexvars, coefs)
    # Add constraints to index.
    consi = MOI.ConstraintIndex{typeof(func), typeof(set)}(num_cons)
    model.constraint_mapping[consi] = num_cons
    return consi
end

function MOI.add_constraint(model::Optimizer,
                            func::MOI.ScalarQuadraticFunction{Float64}, set::LS)
    (model.number_solved >= 1) && throw(AddConstraintError())
    check_inbounds(model, func)
    # We add a constraint in KNITRO.
    num_cons = KN_add_con(model.inner)
    # Add upper bound.
    if isa(set, MOI.LessThan{Float64})
        val = check_value(set.upper)
        KN_set_con_upbnd(model.inner, num_cons, val - func.constant)
    elseif isa(set, MOI.GreaterThan{Float64})
        val = check_value(set.lower)
        KN_set_con_lobnd(model.inner, num_cons, val - func.constant)
    elseif isa(set, MOI.EqualTo{Float64})
        val = check_value(set.value)
        KN_set_con_eqbnd(model.inner, num_cons, val - func.constant)
    end
    # Parse linear structure of constraint.
    indexvars, coefs = canonical_linear_reduction(func)
    KN_add_con_linear_struct(model.inner, num_cons, indexvars, coefs)
    # Parse quadratic term.
    qvar1, qvar2, qcoefs = canonical_quadratic_reduction(func)
    # Take care that Knitro is 0-indexed!
    qvar1 .= qvar1 .- 1
    qvar2 .= qvar2 .- 1

    KN_add_con_quadratic_struct(model.inner, num_cons, qvar1, qvar2, qcoefs)
    # Add constraints to index.
    ci = MOI.ConstraintIndex{typeof(func), typeof(set)}(num_cons)
    model.constraint_mapping[ci] = num_cons
    return ci
end

function MOI.add_constraint(model::Optimizer,
                            func::MOI.VectorAffineFunction, set::MOI.AbstractVectorSet)
    (model.number_solved >= 1) && throw(AddConstraintError())
    # TODO: add check inbounds for VectorAffineFunction.
    previous_col_number = number_constraints(model)
    num_cols = length(func.constants)
    # Add constraints inside KNITRO.
    index_cons = KN_add_cons(model.inner, num_cols)

    # Parse vector affine expression.
    indexcols, indexvars, coefs = canonical_vector_affine_reduction(func)
    # Reformate indexcols with current numbering.
    indexcols .+= previous_col_number

    # Load inside KNITRO.
    KN_add_con_linear_struct(model.inner, indexcols, indexvars, coefs)
    if isa(set, MOI.Nonnegatives)
        KN_set_con_lobnds(model.inner, index_cons, - func.constants)
    elseif isa(set, MOI.Nonpositives)
        KN_set_con_upbnds(model.inner, index_cons, - func.constants)
    elseif isa(set, MOI.Zeros)
        KN_set_con_eqbnds(model.inner, index_cons, - func.constants)
    else
        error("Invalid set $set for VectorAffineFunction constraint")
    end
    # Add constraints to index.
    ci = MOI.ConstraintIndex{typeof(func), typeof(set)}(index_cons[1])
    model.constraint_mapping[ci] = index_cons
    return ci
end

# Add second order cone constraint.
function MOI.add_constraint(model::Optimizer,
                            func::MOI.VectorAffineFunction, set::MOI.SecondOrderCone)
    (model.number_solved >= 1) && throw(AddConstraintError())
    # Add constraints inside KNITRO.
    index_con = KN_add_con(model.inner)

    # Parse vector affine expression.
    indexcoords, indexvars, coefs = canonical_vector_affine_reduction(func)
    constants = func.constants
    # Distinct two parts of secondordercone.
    # First row corresponds to linear part of SOC.
    indlinear = indexcoords .== 0
    indcone = indexcoords .!= 0
    ncoords = length(constants) - 1
    @assert ncoords == set.dimension - 1

    # Load Second Order Conic constraint.
    ## i) linear part
    KN_set_con_upbnd(model.inner, index_con, constants[1])
    KN_add_con_linear_struct(model.inner, index_con,
                             indexvars[indlinear], -coefs[indlinear])

    ## ii) soc part
    index_var_cone = indexvars[indcone]
    nnz = length(index_var_cone)
    index_coord_cone = convert.(Cint, indexcoords[indcone] .- 1)
    coefs_cone = coefs[indcone]
    const_cone = constants[2:end]

    KN_add_con_L2norm(model.inner,
                      index_con, ncoords, nnz,
                      index_coord_cone,
                      index_var_cone,
                      coefs_cone,
                      const_cone)

    # set specific Knitro's params
    KN_set_param(model.inner, KN_PARAM_BAR_CONIC_ENABLE, KN_BAR_CONIC_ENABLE_SOC)
    KN_set_param(model.inner, KN_PARAM_ALGORITHM, KN_ALG_BAR_DIRECT)
    KN_set_param(model.inner, KN_PARAM_BAR_MURULE, KN_BAR_MURULE_FULLMPC)

    # Add constraints to index.
    ci = MOI.ConstraintIndex{typeof(func), typeof(set)}(index_con)
    model.constraint_mapping[ci] = indexvars
    return ci
end

function MOI.add_constraint(model::Optimizer, func::VOV, set::T) where T <: VLS
    (model.number_solved >= 1) && throw(AddConstraintError())
    indv = convert.(Cint, [v.value for v in func.variables] .- 1)
    bnd = zeros(Float64, length(indv))

    if isa(set, MOI.Zeros)
        KN_set_var_fxbnds(model.inner, indv, bnd)
    elseif isa(set, MOI.Nonnegatives)
        KN_set_var_lobnds(model.inner, indv, bnd)
    elseif isa(set, MOI.Nonpositives)
        KN_set_var_upbnds(model.inner, indv, bnd)
    end

    # TODO
    ncons = MOI.get(model, MOI.NumberOfConstraints{VOV, T}())
    ci = MOI.ConstraintIndex{VOV, T}(ncons)
    model.constraint_mapping[ci] = indv
    return ci
end

function MOI.add_constraint(model::Optimizer,
                            func::MOI.VectorOfVariables, set::MOI.SecondOrderCone)
    (model.number_solved >= 1) && throw(AddConstraintError())
    # Add constraints inside KNITRO.
    index_con = KN_add_con(model.inner)
    indv = [v.value - 1 for v in func.variables]

    KN_set_con_upbnd(model.inner, index_con, 0.)
    KN_add_con_linear_struct(model.inner, index_con, indv[1], -1.0)

    indexVars = convert.(Cint, indv[2:end])
    nnz = length(indexVars)
    indexCoords = Cint[i for i in 0:(nnz-1)]
    coefs = ones(Float64, nnz)
    constants = zeros(Float64, nnz)

    KN_add_con_L2norm(model.inner, index_con, nnz, nnz,
                      indexCoords, indexVars, coefs, constants)

    KN_set_param(model.inner, KN_PARAM_BAR_CONIC_ENABLE, KN_BAR_CONIC_ENABLE_SOC)
    KN_set_param(model.inner, KN_PARAM_ALGORITHM, KN_ALG_BAR_DIRECT)
    KN_set_param(model.inner, KN_PARAM_BAR_MURULE, KN_BAR_MURULE_FULLMPC)

    # Add constraints to index.
    ci = MOI.ConstraintIndex{typeof(func), typeof(set)}(index_con)
    model.constraint_mapping[ci] = convert.(Cint, indv)
    return ci
end

# Complementarity constraint

# Complementarity constraints (M x + b  complements x)
# As Knitro supports only complementarity constraints between variables,
# we reformulate it by adding auxiliary variables:
#
# x_aux = Mx + b
# (x_aux complements x)
#
function MOI.add_constraint(model::Optimizer,
                            func::MOI.VectorAffineFunction, set::MOI.Complements)
    (model.number_solved >= 1) && throw(AddConstraintError())

    # Add auxiliary variables `x_aux`.
    x_aux = KN_add_vars(model.inner, set.dimension)

    # Add new constraints in Knitro to formulate x_aux = Mx + b
    n_cons = KN_add_cons(model.inner, set.dimension)
    offset = n_cons[1]

    # Parse VectorAffineFunction defining the complementarity constraints
    indexcoords, indexvars, coefs = canonical_vector_affine_reduction(func)

    # Convert VectorAffineFunction in Knitro format
    x_comp = Cint[]
    jac_cons = Cint[]
    jac_vars = Cint[]
    jac_coefs = Cdouble[]
    for (c, v, f) in zip(indexcoords, indexvars, coefs)
        # By convention, if c is greater than the dimension of
        # the complementarity constraint set, then it specifies
        # the complementarity variable `x`.
        if c >= set.dimension
            push!(x_comp, v)
        else
            # Otherwise, the index corresponds to one of the matrix `M`.
            # Add it to the Jacobian structure of the problem.
            push!(jac_cons, c + offset)
            push!(jac_vars, v)
            push!(jac_coefs, f)
        end
    end
    # Add the index corresponding to `x_aux` in the Jacobian.
    for (icons, ivar) in zip(n_cons, x_aux)
        push!(jac_cons, icons)
        push!(jac_vars, ivar)
        push!(jac_coefs, -1.0)
    end
    # Get constant structure
    q = func.constants[1:set.dimension]

    # Add structure of new constaints to Knitro.
    KN_add_con_linear_struct(model.inner, jac_cons, jac_vars, jac_coefs)
    KN_set_con_eqbnds(model.inner, n_cons, -q)

    # Currently, only complementarity constraints between two variables
    # are supported.
    comp_type = fill(KN_CCTYPE_VARVAR, set.dimension)

    KN_set_compcons(model.inner, comp_type, x_comp, x_aux)
    # Add constraints to index (KNITRO supports only 1 complementarity
    # contraint).
    ci = MOI.ConstraintIndex{typeof(func), typeof(set)}(0)
    return ci
end

# Complementarity constraints (x_1 complements x_2), with x_1 and x_2
# being two variables of the problem.
function MOI.add_constraint(model::Optimizer,
                            func::MOI.VectorOfVariables, set::MOI.Complements)
    (model.number_solved >= 1) && throw(AddConstraintError())
    indv = Cint[v.value - 1 for v in func.variables]
    n_comp = set.dimension
    # Currently, only complementarity constraints between two variables
    # are supported.
    comp_type = fill(KN_CCTYPE_VARVAR, n_comp)

    KN_set_compcons(model.inner, comp_type,
                    indv[1:n_comp],
                    indv[n_comp+1:end])

    # Add constraints to index (KNITRO supports only 1 complementarity
    # contraint).
    ci = MOI.ConstraintIndex{typeof(func), typeof(set)}(0)
    return ci
end

##################################################
## Binary & Integer constraints.

# Define integer and boolean constraints.
function MOI.add_constraint(model::Optimizer,
                            v::MOI.SingleVariable, ::MOI.ZeroOne)
    vi = v.variable
    indv = vi.value - 1
    check_inbounds(model, vi)
    model.number_zeroone_constraints += 1
    # Made the bounds explicit for KNITRO and take care of
    # preexisting MOI's bounds.
    lobnd = 0.
    upbnd = 1.
    KN_set_var_lobnds(model.inner, indv, lobnd)
    KN_set_var_upbnds(model.inner, indv, upbnd)
    KN_set_var_type(model.inner, vi.value - 1, KN_VARTYPE_BINARY)

    return MOI.ConstraintIndex{MOI.SingleVariable, MOI.ZeroOne}(model.number_zeroone_constraints)
end

function MOI.add_constraint(model::Optimizer,
                            v::MOI.SingleVariable, ::MOI.Integer)
    (model.number_solved >= 1) && throw(AddConstraintError())
    vi = v.variable
    check_inbounds(model, vi)
    model.number_integer_constraints += 1
    KN_set_var_type(model.inner, vi.value - 1, KN_VARTYPE_INTEGER)
    return MOI.ConstraintIndex{MOI.SingleVariable, MOI.Integer}(model.number_integer_constraints)
end

##################################################
## Constraint DualStart
function MOI.supports(::Optimizer,
                      ::MOI.ConstraintDualStart,
                      ::MOI.ConstraintIndex{<:SF, <:SS})
    return true
end
function MOI.set(model::Optimizer,
                 ::MOI.ConstraintDualStart,
                 ci::MOI.ConstraintIndex{<:SF,<:SS},
                 value::Union{Real, Nothing})
    if isa(value, Real)
        KN_set_con_dual_init_values(model.inner, ci.value, Cdouble(value))
    else
        # By default, initial value is set to 0.
        KN_set_con_dual_init_values(model.inner, ci.value, Cdouble(0.0))
    end
    return
end

function MOI.supports(::Optimizer,
                      ::MOI.ConstraintDualStart,
                      ::MOI.ConstraintIndex{MOI.SingleVariable, <:LS})
    return true
end
function MOI.set(model::Optimizer,
                 ::MOI.ConstraintDualStart,
                 ci::MOI.ConstraintIndex{MOI.SingleVariable, <:LS},
                 value::Union{Real, Nothing})
    if isa(value, Real)
        KN_set_var_dual_init_values(model.inner, ci.value, Cdouble(value))
    else
        # By default, initial value is set to 0.
        KN_set_var_dual_init_values(model.inner, ci.value, Cdouble(0.0))
    end
    return
end

##################################################
## Constraint naming
# TODO: dry supports with macros.
MOI.supports(::Optimizer, ::MOI.ConstraintName, ::Type{MOI.ConstraintIndex}) = true


text:
cd("../src") do
    run(`make`)
end

ext: jl


text:
const IntVector{w,T} = IntArray{w,T,1}

function IntVector{w,T}(len::Integer, mmap::Bool=false) where {w,T}
    return IntArray{w,T}((len,), mmap)
end

function IntVector{w,T}(mmap::Bool=false) where {w,T}
    return IntArray{w,T}((0,), mmap)
end

function IntVector{w}(vector::AbstractVector{T}) where {w,T}
    return IntArray{w,T,1}(vector)
end

function Base.convert(::Type{IntVector{w}}, vector::AbstractVector{T}) where {w,T}
    return convert(IntArray{w,T,1}, vector)
end

function Base.resize!(vector::IntVector, len::Integer)
    resize!(vector.buffer, len)
    vector.size = (len,)
    return vector
end

function Base.push!(vector::IntVector, x::Integer)
    resize!(vector, length(vector) + 1)
    vector[end] = x
    return vector
end

function Base.pop!(vector::IntVector)
    x = vector[end]
    resize!(vector, length(vector) - 1)
    return x
end

function Base.append!(vec::IntVector, items::AbstractVector)
    len = length(vec)
    resize!(vec, len + length(items))
    for i in 1:lastindex(items)
        vec[len+i] = items[i]
    end
    return vec
end

function Base.reverse!(vec::IntVector, lo::Integer=1, hi::Integer=lastindex(vec))
    return reverse!(vec, Int(lo), Int(hi))
end

function Base.reverse!(vec::IntVector, lo::Int, hi::Int)
    if hi ≤ lo
        return vec
    end
    for i in 0:div(hi - lo, 2)
        vec[lo+i], vec[hi-i] = vec[hi-i], vec[lo+i]
    end
    return vec
end

radixsort(vector::IntVector) = radixsort!(copy(vector))

function radixsort!(vector::IntVector{w}) where {w}
    return radixsort!(vector, 1, length(vector), w)
end

function radixsort!(v::IntVector{w,T}, lo, hi, k) where {w,T}
    @assert 1 ≤ k ≤ w
    if lo > hi
        return v
    end
    i, j = lo, hi
    bit_k = T(1) << (k - 1)
    while i < j
        v_i = v[i]
        if v_i & bit_k == 0
            i += 1
        else
            v[i], v[j] = v[j], v_i
            j -= 1
        end
    end
    @assert i == j
    if k ≥ 2
        if v[i] & bit_k != 0
            i -= 1
        else
            j += 1
        end
        radixsort!(v, lo, i, k - 1)
        radixsort!(v, j, hi, k - 1)
    end
    return v
end


text:
using Revise
using ANOVADDPTest
using Gadfly
using StatsModels
using Statistics
using DataFrames
using Distributions
using Random

# # Simulate sample
# function simulate_sample_normal(rng, N)
#     X = rand(rng, 1:2, N, 2);
#     y = randn(rng, N);
#     for i in 1:N
#         ((X[i, 1] == 2) && (X[i, 2] == 2)) && (y[i] += 1)
#     end
#     return y, X
# end

# Simulate sample
function simulate_sample_poisson(rng, N)
    X = rand(rng, 1:2, N, 2);
    y = rand(rng, Poisson(1.0), N);
    for i in 1:N
        ((X[i, 1] == 2) && (X[i, 2] == 2)) && (y[i] += 1)
    end
    return y, X
end

# Example
N = 1000;
rng = MersenneTwister(1);
y, X = simulate_sample_poisson(rng, N);
fit = anova_bnp_poisson(y, X);
show(group_probs(fit); allrows = true)
show(group_codes(fit); allrows = true)
plot(
    fpost(fit),
    x = :y,
    y = :f,
    color = :group,
    Geom.line,
    Scale.color_discrete()
)


text:
#=

Bit rejection-sampling using discrete
Data Augmentation for Biogeographic Histories.

Ignacio Quintero Mächler

t(-_-t)

May 01 2017

=#

# reprehensible code... but faster
"""
    bit_rejsam!(Y    ::Array{Int64,3},
                idx  ::UnitRange{Int64},
                sf   ::Int64,
                λ1   ::Float64,
                λ0   ::Float64,
                cumts::Array{Float64,1})

Bit rejection-sample a single branch
given the start and end states and a
vector of cumulative δtimes. It assigns to `Y`
in place, avoiding extra memory allocation.
"""
function bit_rejsam!(Y    ::Array{Int64,3},
                     idx  ::UnitRange{Int64},
                     sf   ::Int64,
                     λ1   ::Float64,
                     λ0   ::Float64,
                     cumts::Array{Float64,1})

  idx_end = idx[end]::Int64

  @inbounds begin

    cur_s = Y[idx[1]]::Int64
    cur_t = 0.0
    s     = idx[2]::Int64
    idx_1 = (idx[1] - 1)::Int64

    if cur_s == 0

      while true

        cur_t += rexp(λ1)::Float64
        f      = (idx_1 + idxlessthan(cumts, cur_t))::Int64

        Y[s:f] .= cur_s::Int64

        if f == idx_end
          break
        end

        cur_s = (1 - cur_s)::Int64
        s     = (f + 1)::Int64

        # same but with loss rate
        cur_t += rexp(λ0)::Float64
        f      = (idx_1 + idxlessthan(cumts, cur_t))::Int64

        Y[s:f] .= cur_s::Int64

        if f == idx_end
          break
        end

        cur_s = (1 - cur_s)::Int64
        s     = (f + 1)::Int64

      end

    else

      while true

        cur_t += rexp(λ0)::Float64
        f      = (idx_1 + idxlessthan(cumts, cur_t))::Int64

        Y[s:f] .= cur_s::Int64

        if f == idx_end
          break
        end

        cur_s = (1 - cur_s)::Int64
        s     = (f + 1)::Int64

        # same but with loss rate
        cur_t += rexp(λ1)::Float64
        f      = (idx_1 + idxlessthan(cumts, cur_t))::Int64

        Y[s:f] .= cur_s::Int64

        if f == idx_end
          break
        end

        cur_s = (1 - cur_s)::Int64
        s     = (f + 1)::Int64

      end
    end

  end

  # rejection sampling if end simulation state do not match observed state
  while Y[idx_end] != sf

    cur_s = Y[idx[1]]::Int64
    cur_t = 0.0
    s     = idx[2]::Int64

    if cur_s == 0

      while true

        cur_t += rexp(λ1)::Float64
        f      = (idx_1 + idxlessthan(cumts, cur_t))::Int64

        Y[s:f] .= cur_s::Int64

        if f == idx_end
          break
        end

        cur_s = (1 - cur_s)::Int64
        s     = (f + 1)::Int64

        # same but with loss rate
        cur_t += rexp(λ0)::Float64
        f      = (idx_1 + idxlessthan(cumts, cur_t))::Int64

        Y[s:f] .= cur_s::Int64

        if f == idx_end
          break
        end

        cur_s = (1 - cur_s)::Int64
        s     = (f + 1)::Int64

      end

    else

      while true

        cur_t += rexp(λ0)::Float64
        f      = (idx_1 + idxlessthan(cumts, cur_t))::Int64

        Y[s:f] .= cur_s::Int64

        if f == idx_end
          break
        end

        cur_s = (1 - cur_s)::Int64
        s     = (f + 1)::Int64

        # same but with loss rate
        cur_t += rexp(λ1)::Float64
        f      = idx_1 + idxlessthan(cumts, cur_t)::Int64

        Y[s:f] .= cur_s::Int64

        if f == idx_end
          break
        end

        cur_s = (1 - cur_s)::Int64
        s     = (f + 1)::Int64

      end
    end

  end

end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
# struct SgramData <: JAXTAMData
#     mission::Symbol
#     instrument::Symbol
#     obsid::String
#     bin_time::Real
#     bin_size::Int
#     gti_index::Int
#     gti_start_time::Real
#     amps::Array
#     freqs::Array
# end

function _sgram(lc::JAXTAM.BinnedData, stft_intervals=round(Int, 1024/(lc.bin_time)))
    gti_only_counts = Array{Int64,1}()
    gti_only_times  = Array{Float64,1}()

    for gti in DataFrames.eachrow(lc.gtis)
        gti_start = gti[:start]
        gti_stop  = gti[:stop]

        if gti_stop - gti_start >= 16
            start = findfirst(lc.times .>= gti_start)
            stop  = findfirst(lc.times .>= gti_stop) - 1

            append!(gti_only_counts, lc.counts[start:stop])
            append!(gti_only_times, lc.times[start:stop])
        end
    end

    sgrm = spectrogram(gti_only_counts, stft_intervals, 0; fs=1/lc.bin_time)

    sgram_pseudo_time = 1:size(sgrm.power, 2)
    sgram_freqs  = collect(sgrm.freq)
    sgram_powers = sgrm.power
    sgram_powers[1, :] .= NaN
    sgram_powers[sgram_powers .== -Inf] .= NaN

    heatmap(sgram_freqs, sgram_pseudo_time, sgram_powers', legend=true, size=(1140,400))
end

function _stft(lc::JAXTAM.BinnedData, stft_intervals=round(Int, 1024/(lc.bin_time)))
    gti_only_counts = Array{Int64,1}()
    gti_only_times  = Array{Float64,1}()

    for gti in DataFrames.eachrow(lc.gtis)
        gti_start = gti[:start]
        gti_stop  = gti[:stop]

        if gti_stop - gti_start >= 16
            start = findfirst(lc.times .>= gti_start)
            stop  = findfirst(lc.times .>= gti_stop) - 1

            append!(gti_only_counts, lc.counts[start:stop])
            append!(gti_only_times, lc.times[start:stop])
        end
    end

    dsp_stft      = abs.(stft(gti_only_counts, stft_intervals; fs=1/lc.bin_time))
    #dsp_stft      = dsp_stft.^2
    dsp_stft_time = collect(1:size(dsp_stft, 2)).*(lc.bin_time*stft_intervals/2)
    dsp_stft_freq = linspace(0, 0.5*(1/lc.bin_time), size(dsp_stft, 1))

    dsp_stft[1, :] .= NaN # NaN 0 Hz

    return heatmap(dsp_stft_freq, dsp_stft_time, dsp_stft', legend=true, size=(1140,400))
end


text:
mutable struct MultivariateGammaLikelihood <: MyXMLElement
    el::XMLOrNothing
    data::Parameter
    shapes::Vector{Float64}
    scales::Vector{Float64}
    id::String

    function MultivariateGammaLikelihood(data::Parameter,
                                         shapes::Vector{Float64},
                                         scales::Vector{Float64},
                                         id::String)
        @assert length(data.val) == length(shapes) == length(scales)
        return new(nothing, data, shapes, scales, id)
    end
end

function make_xml(mgl::MultivariateGammaLikelihood)
    make_xml(mgl.data)

    el = new_element(bn.MULTIVARIATE_GAMMA_LIKELIHOOD)
    set_id!(el, mgl.id)

    data_el = new_child(el, bn.DATA)
    add_ref_el(data_el, mgl.data.el)

    shape_el = new_child(el, bn.SHAPE)
    add_parameter(shape_el, value = mgl.shapes, lower = 0.0)

    scale_el = new_child(el, bn.SCALE)
    add_parameter(scale_el, value = mgl.scales, lower = 0.0)

    mgl.el = el
    return el
end

function get_precision_prior(mgl::MultivariateGammaLikelihood)
    make_xml(mgl)
    return mgl.el
end

function set_shrinkage_mults!(mgl::MultivariateGammaLikelihood;
                            shapes::Vector{Float64} = Float64[],
                            scales::Vector{Float64} = Float64[])

    n = length(mgl.shapes)
    set_values(shapes, mgl.shapes)
    set_values(scales, mgl.scales)
end

function set_values(src::Vector{Float64}, dest::Vector{Float64})

    if length(src) > 0
        if length(src) == length(dest) - 1
            dest[2:end] .= src
        else
            dest .= src
        end
    end
end


text:
# Convert piece-wise constant p (design region) to pvec (whole domain)
function p_vec(p,P,tags,design_tag)
    pvec = zeros(num_free_dofs(P))
    pi = 0
    @assert length(tags)==num_free_dofs(P)
    for i=1:length(tags)
        if tags[i] == design_tag
            pi += 1
            pvec[i] = p[pi]
        end
    end
    pvec
end

# Extract the design region part from a whole vector
function extract_design(pvec,np,tags,design_tag)
    p_d = zeros(eltype(pvec),np)
    pi = 0
    @assert length(pvec)==length(tags)
    for i=1:length(tags)
        if tags[i] == design_tag
            pi += 1
            p_d[pi] = pvec[i]
        end
    end
    @assert np==pi
    p_d
end

# Gaussian Distribution function with center x0
function GaussianD(x,x0::AbstractArray,δ::AbstractArray)
    n=length(x)
    @assert (n==length(x0))&&(n==length(δ))
    δn = 1.0
    x_δ = 0.0
    for i=1:n
        δn *= √(2π)*δ[i]
        x_δ += ((x[i]-x0[i])/δ[i])^2
    end
    1.0/δn*exp(-x_δ/2.0)
end

# Gaussian Distribution function with center x0
function GaussianY(x,x0::AbstractArray,δ::AbstractArray)
    n=length(x)
    @assert (n==length(x0))&&(n==length(δ))
    δn = √(2π)*δ[2]
    x_δ = ((x[2]-x0[2])/δ[2])^2
    1.0/δn*exp(-x_δ/2.0)
end

max_stars_repo_name: Tiamat-Tech/hpinn

max_issues_repo_name: louhz/hpinn

max_forks_repo_name: louhz/hpinn


text:
using Distributions
using Base.Test
import Distributions: BinomialGeomSampler, BinomialTPESampler, BinomialPolySampler, BinomialAliasSampler

function test_binomsampler(s, params::(Int, Float64), ns::Int, tol::Float64)
    # s: the sampler to be tested
    # params: the binomial parameters, in the form of (n, p)
    # n: the number of samples to be generated
    # tol: the tolerance of deviation

    n, p = params
    if p == 0.0 || n == 0
        for i=1:ns
            @assert rand(s) == 0
        end
    elseif p == 1.0
        for i=1:ns
            @assert rand(s) == n
        end
    else
        cnts = zeros(Int, n+1)
        for i=1:ns
            x = rand(s)
            @assert 0 <= x <= n
            cnts[x+1] += 1
        end
        pv = Distributions.binompvec(n, p)
        pr = cnts ./ ns
        @test_approx_eq_eps pr pv tol
    end
end

# test cases

for params in [(0, 0.4), (0, 0.6), (5, 0.0), (5, 1.0), (1, 0.2), (1, 0.8), (3, 0.4), (4, 0.6)]
    test_binomsampler(BinomialGeomSampler(params...), params, 10^5, 0.015)
end

for params in [(40, 0.5), (100, 0.4), (300, 0.6)]
    test_binomsampler(BinomialTPESampler(params...), params, 10^6, 0.015)
end

for params in [(0, 0.4), (0, 0.6), (5, 0.0), (5, 1.0),
               (1, 0.2), (1, 0.8), (3, 0.4), (4, 0.6),
               (40, 0.5), (100, 0.4), (300, 0.6)]
    test_binomsampler(BinomialPolySampler(params...), params, 10^6, 0.015)
end

for params in [(0, 0.4), (0, 0.6), (5, 0.0), (5, 1.0),
               (1, 0.2), (1, 0.8), (3, 0.4), (4, 0.6),
               (40, 0.5), (100, 0.4), (300, 0.6)]
    test_binomsampler(BinomialAliasSampler(params...), params, 10^6, 0.015)
end


text:
# Parameter helper functions
const gas_molpL = 22.4

watts_to_light_mol(watts) = watts * 4.57e-6
light_mol_to_watts(light_mol) = light_mol / 4.57e-6
water_content_to_mols_per_litre(wc) = wc * 55.5 # L/L of water to mol/L
fraction_per_litre_gas_to_mols(frac) = frac / 22.4

"""
    CarbonVars(J_L_F, X_C, X_O, soilwaterpotential)

Variables for carbon assimilation

$(FIELDDOCTABLE)
"""
@flattenable @udefault_kw @units @bounds @description mutable struct CarbonVars{MoMS,MoL,KPA}
    # Field                  | Flatn | Default i                 | Unit          | Bounds                             | Description
    J_L_F::MoMS              | false | watts_to_light_mol(800.0) | mol*m^-2*s^-1 | watts_to_light_mol.((0.0, 2000.0)) | "Flux of useful photons"
    X_C::MoL                 | false | 400.0*1e-6                | mol*L^-1      | (200.0, 700.0) .* 1e-6             | "Carbon dioxide concentration in air"
    X_O::MoL                 | false | 0.21 * gas_molpL          | mol*L^-1      | (0.1, 0.3) .* gas_molpL            | "Oxygen concentration in air"
    soilwaterpotential::KPA  | false | -100.0                    | kPa           | (0.0, -10000.0)                    | "Soil water potential"
end

"""
    NitrogenVars(soilwaterpotential, soilwaterconent, X_NH, X_NO, X_H)

Variables for nitrogen assimilation.

$(FIELDDOCTABLE)
"""
@flattenable @udefault_kw @units @bounds @description mutable struct NitrogenVars{F,KPA,MoL}
    soilwaterpotential::KPA  | false | -100.0 | kPa      | (0.0, -10000.0) | "Soil water potential"
    soilwaterconent::F       | false | -100.0 | _        | (0.0, -10000.0) | _
    X_NH::MoL                | false | 0.005  | mol*L^-1 | (0.0, 0.1)      | "Concentration of ammonia"
    X_NO::MoL                | false | 0.01   | mol*L^-1 | (0.0, 0.1)      | "Concentration of nitrate see e.g. (_@crawford1998molecular)"
    X_H::MoL                 | false | 10.0   | mol*L^-1 | (0.0, 20.0)     | _
end

" Abstract supertype of all Assimilation components"
abstract type AbstractAssim end

" Abstract supertype of all Carbon assimilation types"
abstract type AbstractCAssim <: AbstractAssim end

" Abstract supertype of all Nitrogen assimilation types"
abstract type AbstractNAssim <: AbstractAssim end

" Abstract supertype of Ammonia/Nitrate assimilation types"
abstract type AbstractNH4_NO3Assim <: AbstractNAssim end

"""
    ConstantCAssim(c_uptake)

C is assimilated at a constant rate, without control from the environment

$(FIELDDOCTABLE)
"""
@columns struct ConstantCAssim{μMoMS} <: AbstractCAssim
    # Field         | Def | Unit              | Bounds      | Log | Description
    c_uptake::μMoMS | 0.1 | μmol*mol^-1*s^-1  | (0.0, 10.0) | _   | "Constant rate of C uptake"
end

"""
    ConstantNAssim(n_uptake)

N is assimilated at a constant rate, without control from the environment

$(FIELDDOCTABLE)
"""
@columns struct ConstantNAssim{μMoMS} <: AbstractNAssim
    n_uptake::μMoMS | 0.1  | μmol*mol^-1*s^-1 | (0.0, 0.5)  | _   | "Constant rate of N uptake"
end

@mix @columns struct SLA{MG}
    SLA::MG         | 24.0 | m^2*kg^-1        | (5.0, 30.0) | _   | "Specific leaf Area. Ferns: 17.4, Forbs: 26.2, Graminoids: 24.0, Shrubs: 9.10, Trees: 8.30"
end

# Mixin paramters for multiple similar formulations
@mix @columns @SLA struct MixinKooijmanPhoto{μMoMoS,MoL,μMoMS,MoMS,V}
    #Field              | Default           | Unit             | Bounds           | Log  | Description
    vars::V             | CarbonVars()      | _                | _                | _    | _
    k_C_binding::μMoMoS | 10000.0           | μmol*mol^-1*s^-1 | (1e-5, 2000.0)   | true | "Scaling rate for carbon dioxide"
    k_O_binding::μMoMoS | 10000.0           | μmol*mol^-1*s^-1 | (1e-5, 2000.0)   | true | "Scaling rate for oxygen"
    K_C::MoL            | 50*1e-6/gas_molpL | mol*L^-1         | (1e-7, 100.0)    | true | "Half-saturation concentration of carbon dioxide"
    K_O::MoL            | 0.0021/gas_molpL  | mol*L^-1         | (1e-7, 100.0)    | true | "Half-saturation concentration of oxygen"
    J_L_K::MoMS         | 2000.0            | mol*m^-2*s^-1    | (1e-3, 100000.0) | true | "Half-saturation flux of useful photons"
    j_L_Amax::μMoMS     | 100.01            | μmol*m^-2*s^-1   | (1e-4, 10000.0)  | true | "Maximum specific uptake of useful photons"
    j_C_Amax::μMoMS     | 20.0              | μmol*m^-2*s^-1   | (5.0,  100.0)    | true | "Maximum specific uptake of carbon dioxide"
    j_O_Amax::μMoMS     | 0.1               | μmol*m^-2*s^-1   | (0.01,  50.0)    | true | "Maximum specific uptake of oxygen"
end

abstract type AbstractKooijmanPhoto <: AbstractCAssim end

"""
    KooijmanSLAPhotosynthesis(vars, k_C_binding, k_O_binding, K_C, K_O, J_L_K, j_L_Amax, j_C_Amax, j_O_Amax)

Parameters for simple photosynthesis module. With specific leaf area to convert area to mass

$(FIELDDOCTABLE)
"""
@MixinKooijmanPhoto struct KooijmanSLAPhotosynthesis{} <: AbstractKooijmanPhoto end

# @columns struct WaterPotentialCutoff{KPA}
    # cutoff::KPA    | -250.0  | kPa  | (0.0, -500.0) | _ | "Max spec uptake of oxygen"
# end

"""
Parameters for Ammonia/Nitrate assimilation

$(FIELDDOCTABLE)
"""
@columns struct KooijmanNH4_NO3Assim{μMoMS,F,MoMo,MoL,V} <: AbstractNH4_NO3Assim
    vars::V          | NitrogenVars() | _         | _             | _ | _
    j_NH_Amax::μMoMS | 50.0    | μmol*mol^-1*s^-1 | (0.1, 1000.0) | _ | "Max spec uptake of ammonia"
    j_NO_Amax::μMoMS | 50.0    | μmol*mol^-1*s^-1 | (0.1, 1000.0) | _ | "Max spec uptake of nitrate"
    ρNO::F           | 0.7     | _                | (1e-4, 1.0)   | _ | "Weights preference for nitrate relative to ammonia." # 1 or less but why?
    y_E_CH_NH::MoMo  | 1.25    | mol*mol^-1       | (1e-4, 4.0)   | _ | "From roots C-reserve to reserve using ammonia"
    K_NH::MoL        | 0.01    | mol*L^-1         | (1e-4, 10.0)  | _ | "Half-saturation concentration of ammonia"
    K_NO::MoL        | 0.01    | mol*L^-1         | (1e-4, 10.0)  | _ | "Half-saturation concentration of nitrate"
    K_H::MoL         | 10.0    | mol*L^-1         | (5.0, 20.0)   | _ | "Half-saturation concentration of water"
end

"""
Parameters for Nitrogen assimilation

$(FIELDDOCTABLE)
"""
@columns struct NAssim{μMoS,MoL,V} <: AbstractNAssim
    vars::V          | NitrogenVars() | _         | _             | _ | _
    j_N_Amax::μMoS   | 50.0    | μmol*mol^-1*s^-1 | (0.1, 1000.0) | _ | "Max spec uptake of ammonia"
    K_N::MoL         | 0.01    | mol*L^-1         | (1e-4, 1.0)   | _ | "Half-saturation concentration of nitrate"
    K_H::MoL         | 10.0    | mol*L^-1         | (1e-2, 100.0) | _ | "Half-saturation concentration of water"
end

"""
    assimilation!(o::AbstractOrgan, u)

Runs assimilation methods, depending on formulation and state `u`.
"""
function assimilation! end
assimilation!(organs::Tuple, u) = map(assimilation!, organs, u)
assimilation!(o, u) = begin
    isgerminated(o, u) && assimilation!(has_reserves(o), assimilation_pars(o), o, u)
    nothing
end
assimilation!(::Nothing, x, o, u, env) = nothing
assimilation!(::Any, f::AbstractCAssim, o, u) = begin
    c = photosynthesis(f, o, u)
    flux(o)[:C,:asi] = max(c, zero(c)) * u[:V] * scaling(o)
end
assimilation!(p::HasCN, f::AbstractNH4_NO3Assim, o, u) = begin
    J = flux(o)
    J_N_ass, J_NO_ass, J_NH_ass = nitrogen_uptake(f, o, u) .* u[:V] * scaling(o)

    θNH = J_NH_ass/J_N_ass                          # Fraction of ammonia in arriving N-flux
    θNO = 1 - θNH                                   # Fraction of nitrate in arriving N-flux
    y_E_EC = θNH * f.y_E_CH_NH + θNO * shared(o).y_E_EC  # Yield coefficient from C-reserve to reserve

    C_tra = J[:C,:tra]

    # Merge rejected C from shoot and uptaken N into reserves
    J[:C,:tra], J[:N,:asi], J[:E,:asi] = stoich_merge(C_tra, J_N_ass, y_E_EC, 1/n_N_E(o))
    # stoich_merge_losses(c_tra, J_N_ass, J[:C,:tra], J[:N,:asi], J[:E,:asi], 1, 1, n_N_E(o))

    # Unused NH₄ remainder is lost so we recalculate N assimilation for NO₃ only
    J[:N,:asi] = (J_NO_ass - θNO * n_N_E(o) * J[:E,:asi]) * 1/n_N_EN(o)
end
"""
    assimilation!(f::AbstractNAssim, o, u)

Runs nitrogen uptake, and combines N with translocated C.
"""
assimilation!(::Any, f::AbstractNAssim, o, u) = begin
    n = nitrogen_uptake(f, o, u)
    flux(o)[:N,:asi] = max(n, zero(n)) * u[:V] * scaling(o)
end

"""
    photosynthesis(f::AbstractCAssim, o, u)

Return the assimilated C for `Organ` `o` with state variables `u`.
"""
function photosynthesis end

photosynthesis(f::ConstantCAssim, o, u) = f.c_uptake

photosynthesis(f::KooijmanSLAPhotosynthesis, o, u) = begin
    v = vars(o); va = f.vars
    mass_area_coef = w_V(o) * f.SLA

    # Photon flux is not temperature dependent, but O and C flux is.
    # See "Stylized facts in microalgal growth" Lorena, Marques, Kooijman, & Sousa.
    j1_l = half_saturation(f.j_L_Amax, f.J_L_K, va.J_L_F) * mass_area_coef / 2
    j1_c = half_saturation(f.j_C_Amax, f.K_C, va.X_C) * mass_area_coef * tempcorrection(o)
    j1_o = half_saturation(f.j_O_Amax, f.K_O, va.X_O) * mass_area_coef * tempcorrection(o)

    # photorespiration.
    bound_o = j1_o/f.k_O_binding # mol/mol
    bound_c = j1_c/f.k_C_binding # mol/mol

    # C flux
    j_c_intake = (j1_c - j1_o)
    j1_co = j1_c + j1_o
    co_l = j1_co/j1_l - j1_co/(j1_l + j1_co)

    j_c_intake / (1 + bound_c + bound_o + co_l)
end

"""
    nitrogen_uptake(f::ConstantNAssim, o, u)

Returns constant nitrogen assimilation.
"""
nitrogen_uptake(f::ConstantNAssim, o, u) = f.n_uptake * tempcorrection(o)

"""
    nitrogen_uptake(f::KooijmanNH4_NO3Assim, o, u)

Returns total nitrogen, nitrate and ammonia assimilated in mols per time.
"""
function nitrogen_uptake(f::KooijmanNH4_NO3Assim, o, u)
    va = f.vars

    K1_NH = half_saturation(f.K_NH, f.K_H, va.X_H) # Ammonia saturation. va.X_H was multiplied by ox.scaling. But that makes no sense.
    K1_NO = half_saturation(f.K_NO, f.K_H, va.X_H) # Nitrate saturation
    J1_NH_ass = half_saturation(f.j_NH_Amax, K1_NH, va.X_NH) * tempcorrection(o) # Arriving ammonia mols.mol⁻¹.s⁻¹
    J_NO_ass = half_saturation(f.j_NO_Amax, K1_NO, va.X_NO) * tempcorrection(o) # Arriving nitrate mols.mol⁻¹.s⁻¹

    J_N_ass = J1_NH_ass + f.ρNO * J_NO_ass # Total arriving N flux
    return (J_N_ass, J_NO_ass, J1_NH_ass)
end

"""
    nitrogen_uptake(f::NAssim, o, u)

Returns nitrogen assimilated in mols per time.
"""
function nitrogen_uptake(f::NAssim, o, u)
    va = f.vars
    # Ammonia proportion in soil water
    K1_N = half_saturation(f.K_N, f.K_H, va.X_H)
    # Arriving ammonia in mol mol^-1 s^-1
    half_saturation(f.j_N_Amax, K1_N, va.X_NO) * tempcorrection(o)
end


text:
function acquire(num, sem)
    sleep(rand())
    println("Task $num waiting for semaphore")
    lock(sem)
    println("Task $num has acquired semaphore")
    sleep(rand())
    unlock(sem)
end

function runsem(numtasks)
    println("Sleeping and running $numtasks tasks.")
    sem = Base.Threads.RecursiveSpinLock()
    @sync(
    for i in 1:numtasks
        @async acquire(i, sem)
    end)
    println("Done.")
end

runsem(4)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using JLD
using Distances

if !isdefined(:loaded) || !loaded
    println("Loading...")
    jldopen("85.jld") do file
        global X, vocab, vocab2idx, idx2vocab
        X = read(file, "X")
        vocab = read(file, "vocab")
        vocab2idx = read(file, "vocab2idx")
        idx2vocab = read(file, "idx2vocab")
        global loaded = true
    end
end

# note: depends some globals
function similar_words(word; v = nothing, K=10)
    if v === nothing
        v = X[vocab2idx[word],:]
    end
    distances = zeros(length(vocab))
    for i in 1:length(vocab)
        distances[i] = cosine_dist(v, X[i,:])
    end

    p = sortperm(distances)
    words = []
    for i in 1:K
        push!(words, (idx2vocab[p[i]], distances[p[i]]))
    end

    words
end

v1 = X[vocab2idx["spain"],:]
v2 = X[vocab2idx["madrid"],:]
v3 = X[vocab2idx["athens"],:]

v = v1 - v2 + v3

words = similar_words("Dummy", v = v, K=11)

for (idx, (word, score)) in enumerate(words)
    println("$idx: ($word, $score)")
end

#=
julia> @time include("./89.jl")
1: (Spain, 0.1298535311445519)
2: (Sweden, 0.15812807697715514)
3: (Italy, 0.15933862165878132)
4: (Germany, 0.18661154954426418)
5: (Austria, 0.19796824010990155)
6: (France, 0.199986843038379)
7: (Netherlands, 0.20552766319824411)
8: (Belgium, 0.22052314063678613)
9: (Denmark, 0.2216631629606418)
10: (Télévisions, 0.22247732061479797)
11: (Hungary, 0.2550573826225834)
  0.806980 seconds (2.53 M allocations: 1.045 GiB, 9.11% gc time)
=#

# ***90***
#=
julia> @time include("./89.jl")
1: (Spain, 0.3457379922084114)
2: (Greece, 0.4379238679036045)
3: (Sweden, 0.44022662117009537)
4: (Athens, 0.4428090455580942)
5: (Norway, 0.46073073622093297)
6: (Poland, 0.48033036303527465)
7: (Denmark, 0.4824352945520548)
8: (Netherlands, 0.4904862473474464)
9: (Turkey, 0.49091823018130454)
10: (Germany, 0.5010240875959475)
11: (Italy, 0.5093253306302299)
=#


text:
using FEMBase, Test

# From the beginning of a project we had a clear concept in our mind: "everything
# is a field". That is, everything can vary temporally and spatially. We think
# that constant is just a special case of field which does not vary in temporal
# nor spatial direction. Fields can vary in spatial direction, i.e. can be either
# constant or variable, and in temporal direction, i.e. can be time variant or
# time invariant. From this pondering we can think that there exists four kind of
# (discrete) fields:

# - discrete, constant, time invariant (DCTI)
# - discrete, variable, time invariant (DVTI)
# - discrete, constant, time variant (DCTV)
# - discrete, variable, time variant (DVTV)

# Discrete, in this context, means that field is defined in point-wise in
# $1 \ldots n$ locations, from where it is then interpolated to whole domain
# using some interpolation polynomials, i.e.
# ```math
# u(\xi, t) = \sum_{i} u_i[t] N_{i}(\xi,t),
# ```math
# where
#     $N_{i}(\xi, t)$
# is the basis function or interpolation polymial corresponding to $i$^{th}
# discrete value and
#     $u_{i}$
# is the discrete value.

# Then we have continuous fields, which are defined in whole domain, or at least
# not point-wise. By following the already used abbreviations, we have four more
# fields:

# - continuous, constant, time invariant (CCTI)
# - continuous, variable, time invariant (CVTI)
# - continuous, constant, time variant (DCTV)
# - continuous, variable, time variant (CVTV)

# Continuous, again in this context, does not mean that field has to be defined
# everywhere. It's enough that it's defined in function of spatial and/or temporal
# coordinates, i.e. we have $u \equiv u(\xi, t)$, without a some spesific basis
# needed to interpolate from discrete values.

# Field itself can be in principle anything. However, usually either scalar,
# vector or tensor (matrix). Time does not to have be real, it can be for example
# angle of some rotating machine or even complex value.

# From these starting points, we assume that the mentioned field system can
# describe all imaginable situations.

# ## Creating new fields

# For discrete fields that are varying in spatial direction, value for each
# discrete point is defined using NTuple. The order of points is implicitly
# assumed to be same than node ordering in ABAQUS. That is, first corner nodes
# in anti-clockwise direction and after that middle nodes.

# For example, `(1, 2, 3, 4)` is a scalar field having length of 4 and
# `([1,2],[2,3],[3,4],[4,5])` is a vector field having length of 4.

# For fields that are varying in temporal direction, `time => value` syntax is
# used. The first item in pair is time (or similar) and second item is value
# assigned to that time. For example, `0.0 => 1.0` is a time-dependent scalar
# field having value 1.0 at time 0.0.

# ## Dicrete, constant, time invariant field (DCTI)

# The most simple field is a field that is constant in both time and spatial
# direction. Discrete, constant, time invariant field. For example, youngs
# modulus could be this kind of field.

a = DCTI(1)

# Accessing data is done using `interpolate`. In FEM codes, we try to hide the
# actual type of the field, so for example interpolating constant field works,
# but the result is quite unsuprising.

@test interpolate(a, 0.0) == 1

# Field value value can be updated with `update!` function:

update!(a, 2)
@test a == 2

# Constant field of course doesn't have to be scalar field. It can be e.g.
# vector field. I use here packate Tensors.jl because of its excellent
# performance and other features, but normal `Vector` would work just fine
# also:

using Tensors

b = DCTI(Vec(1, 2))

# Interpolation, again, returns just the original data:

@test interpolate(b, 0.0) == [1, 2]

# Updating field is done using `update!`-function:

update!(b, Vec(2, 3))
@test interpolate(b, 0.0) == [2, 3]

# Constant tensor field:

c = DCTI(Tensor{2,2}((1.0, 2.0, 3.0, 4.0)))

# Data can be accessed also using `getindex`. Also things like `length` and
# `size` are defined.

@test interpolate(c, 0.0) == [1 3; 2 4]
@test c[1] == [1 3; 2 4]
@test length(c) == 4
@test size(c) == (2, 2)

# For now everything might look like extra complexity, but later on we see how
# to combine field with some basis functions in order to interpolate in element
# domain. Another nice feature is that we can interpolate fields in time. In this
# particular case of time invariant fields it of course doesn't give anything
# extra.

# ## Dicrete, variable, time invariant fields (DVTI)

@testset "DVTI field" begin

    # scalar field
    a = DVTI((1, 2))
    @test a[1] == 1
    @test a[2] == 2
    @test interpolate(a, 0.0) == (1, 2)
    update!(a, (2, 3))
    @test a == (2, 3)
    @test (2, 3) == a

    # vector field
    b = DVTI(([1, 2], [2, 3]))
    @test b[1] == [1, 2]
    @test b[2] == [2, 3]
    @test interpolate(b, 0.0) == ([1, 2], [2, 3])
    update!(b, ([2, 3], [4, 5]))
    @test b == ([2, 3], [4, 5])

    # tensor field
    c = DVTI(([1 2; 3 4], [2 3; 4 5]))
    @test c[1] == [1 2; 3 4]
    @test c[2] == [2 3; 4 5]
    @test interpolate(c, 0.0) == ([1 2; 3 4], [2 3; 4 5])
    update!(c, ([2 3; 4 5], [5 6; 7 8]))
    @test c == ([2 3; 4 5], [5 6; 7 8])

    d = DVTI(2, 3)
    @test a == d
end

@testset "DCTV field" begin

    # scalar field
    a = DCTV(0.0 => 0.0, 1.0 => 1.0)
    @test isapprox(interpolate(a, -1.0), 0.0)
    @test isapprox(interpolate(a, 0.0), 0.0)
    @test isapprox(interpolate(a, 0.5), 0.5)
    @test isapprox(interpolate(a, 1.0), 1.0)
    update!(a, 1.0 => 2.0)
    @test isapprox(interpolate(a, 0.5), 1.0)
    update!(a, 2.0 => 1.0)
    @test isapprox(interpolate(a, 1.5), 1.5)

    # vector field
    b = DCTV(0.0 => [1.0, 2.0], 1.0 => [2.0, 3.0])
    @test isapprox(interpolate(b, 0.5), [1.5, 2.5])

    # tensor field
    c = DCTV(0.0 => [1.0 2.0; 3.0 4.0], 1.0 => [2.0 3.0; 4.0 5.0])
    @test isapprox(interpolate(c, 0.5), [1.5 2.5; 3.5 4.5])
end

@testset "DVTV field" begin
    # scalar field
    a = DVTV(0.0 => (0.0, 1.0), 1.0 => (1.0, 0.0))
    update!(a, 2.0 => (2.0, 0.0))
    r = interpolate(a, 0.5)
    @test isapprox(r[1], 0.5)
    @test isapprox(r[2], 0.5)
    update!(a, 2.0 => (4.0, 0.0))
end

@testset "CVTV field" begin
    f = CVTV((xi, t) -> xi[1] * xi[2] * t)
    @test isapprox(f([1.0, 2.0], 3.0), 6.0)
end

@testset "Dictionary fields" begin
    X = Dict(1 => [0.0, 0.0], 1000 => [1.0, 0.0], 100000 => [1.0, 1.0])
    G = DVTId(X)
    @test isapprox(G[1], X[1])
    @test isapprox(G[1000], X[1000])
    @test isapprox(G[100000], X[100000])
    Y = Dict(1 => [2.0, 2.0], 1000 => [3.0, 2.0], 100000 => [3.0, 3.0])
    F = DVTVd(0.0 => X, 1.0 => Y)
    @test isapprox(interpolate(F, 0.5)[100000], [2.0, 2.0])
end

@testset "update dictionary field" begin
    f1 = Dict(1 => 1.0, 2 => 2.0, 3 => 3.0)
    f2 = Dict(1 => 2.0, 2 => 3.0, 3 => 4.0)
    fld = DVTVd(0.0 => f1)
    update!(fld, 1.0 => f2)
    @test isapprox(interpolate(fld, 0.5)[1], 1.5)
    update!(fld, 1.0 => f1)
    @test isapprox(interpolate(fld, 0.5)[1], 1.0)
end

@testset "use of common constructor field" begin
    @test isa(field(1.0), DCTI)
    @test isa(field(1.0 => 1.0), DCTV)
    @test isa(field((1.0, 2.0)), DVTI)
    @test isa(field(1, 2), DVTI)
    @test isa(field(1.0 => (1.0, 2.0)), DVTV)
    @test isa(field((xi, t) -> xi[1] * t), CVTV)
    @test isa(field(1 => [1.0, 2.0], 10 => [2.0, 3.0]), DVTId)
    @test isa(field(0.0 => (1 => 1.0, 10 => 2.0), 1.0 => (1 => 2.0, 10 => 3.0)), DVTVd)
    X = Dict(1 => [0.0, 0.0], 2 => [1.0, 0.0])
    X1 = field(X)
    X2 = field(0.0 => X)
    @test isa(X1, DVTId)
    @test isa(X2, DVTVd)
end

@testset "general interpolation" begin
    a = [1, 2, 3]
    b = (2, 3, 4)
    @test interpolate(a, b) == 2 + 6 + 12
    a = (1, 2)
    b = (2, 3, 4)
    @test interpolate(a, b) == 2 + 6
    @test_throws AssertionError interpolate(b, a)
end


text:
using PlotlyJS

"""
`boxplot( val::Vector, cat::Vector; ...)`

Description
===========

Generates a boxplot using numeric data provided by an Array
`val` and categorical (numeric or string) data provided by
`cat`

Arguments
=========

- **`boxpoints`** : ["outliers"], "all", "suspectedoutliers"
- **`boxmean`** : ["True"], "sd", "False"
- **`hovermode`** : ["closest"], "x", "y", "False"
- **`jitter`** : [0.3] Between 0-1
- **`pointpos`** : Numeric
- **`orientation`** : ["v"], "h"
- **`maxcategories`** : threshold to ensure the cardinality of categories isn't too high

"""
function boxplot{N<:Number}( val::Vector{N};
                             cat::Vector=[],
                             title="", height=600, width=600,
                             xaxis_title="", yaxis_title="",
                             boxpoints="outliers", boxmean="True",
                             hovermode="closest", jitter=0.3, pointpos=0,
                             orientation="v", maxcategories=16)

    # If cat is not specified then create one
    if cat == []
        cat = ["0" for i in 1:length(val)]
    end

    # Ensure that cat and val are of equal length
    @assert length(cat) == length(val) "cat and val should be equal in length"

    # Check for too many categories; more than 16 are too many
    @assert length(unique(cat)) <= maxcategories "Too many categories, try a different plot"

    # Check if x is a subtype of Number
    @assert eltype(val) <: Number "x should be of Numeric Type"

    data = Array{PlotlyJS.GenericTrace{Dict{Symbol,Any}},1}([])

    # TODO: Mechanism to specify the order of boxes
    #       For now just sort the categories to have an order
    boxids = sort(unique(cat))

    for id in boxids
        sub = val[cat .== id]
        name=string(id)

        # TODO: A clean method to handle horizontal / vertical box plots
        #       by replacing "y=sub" with "x=sub" will result in a horizontal plot

        if orientation=="h"
            trace = box(;x=sub, name=name, boxpoints=boxpoints,
                        jitter=jitter, pointpos=pointpos,
                        boxmean=boxmean)
        elseif orientation=="v"
            trace = box(;y=sub, name=name, boxpoints=boxpoints,
                        jitter=jitter, pointpos=pointpos,
                        boxmean=boxmean)
        end

        push!(data, trace)
    end

    layout = Layout(;title=title,
                    xaxis_title=xaxis_title, yaxis_title=yaxis_title,
                    width=width, height=height, hovermode=hovermode)

    p = Plot(data, layout)

end


text:
using SafeTestsets

@safetestset "Utility Functions                                                               " begin include("utils_tests.jl") end
@safetestset "Abstract Equation                                                               " begin include("geometric_equations_tests.jl") end
@safetestset "Ordinary Differential Equations                                                 " begin include("ordinary_differential_equations_tests.jl") end
@safetestset "Differential Algebraic Equations                                                " begin include("differential_algebraic_equations_tests.jl") end
@safetestset "Stochastic Differential Equations                                               " begin include("stochastic_differential_equations_tests.jl") end
@safetestset "Abstract Problem                                                                " begin include("geometric_problem_tests.jl") end
@safetestset "Geometric Problem                                                               " begin include("geometric_problem_tests.jl") end


text:
function Vector3D(arg0::Vector{jdouble})
    return Vector3D((Vector{jdouble},), arg0)
end

function Vector3D(arg0::jdouble, arg1::Vector3D)
    return Vector3D((jdouble, Vector3D), arg0, arg1)
end

function Vector3D(arg0::jdouble, arg1::Vector3D, arg2::jdouble, arg3::Vector3D)
    return Vector3D((jdouble, Vector3D, jdouble, Vector3D), arg0, arg1, arg2, arg3)
end

function Vector3D(arg0::jdouble, arg1::Vector3D, arg2::jdouble, arg3::Vector3D, arg4::jdouble, arg5::Vector3D)
    return Vector3D((jdouble, Vector3D, jdouble, Vector3D, jdouble, Vector3D), arg0, arg1, arg2, arg3, arg4, arg5)
end

function Vector3D(arg0::jdouble, arg1::Vector3D, arg2::jdouble, arg3::Vector3D, arg4::jdouble, arg5::Vector3D, arg6::jdouble, arg7::Vector3D)
    return Vector3D((jdouble, Vector3D, jdouble, Vector3D, jdouble, Vector3D, jdouble, Vector3D), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function Vector3D(arg0::jdouble, arg1::jdouble)
    return Vector3D((jdouble, jdouble), arg0, arg1)
end

function Vector3D(arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return Vector3D((jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function add(obj::Vector3D, arg0::JVector)
    return jcall(obj, "add", Vector3D, (JVector,), arg0)
end

function add(obj::Vector3D, arg0::jdouble, arg1::JVector)
    return jcall(obj, "add", Vector3D, (jdouble, JVector), arg0, arg1)
end

function angle(::Type{Vector3D}, arg0::Vector3D, arg1::Vector3D)
    return jcall(Vector3D, "angle", jdouble, (Vector3D, Vector3D), arg0, arg1)
end

function cross_product(::Type{Vector3D}, arg0::Vector3D, arg1::Vector3D)
    return jcall(Vector3D, "crossProduct", Vector3D, (Vector3D, Vector3D), arg0, arg1)
end

function cross_product(obj::Vector3D, arg0::JVector)
    return jcall(obj, "crossProduct", Vector3D, (JVector,), arg0)
end

function distance(::Type{Vector3D}, arg0::Vector3D, arg1::Vector3D)
    return jcall(Vector3D, "distance", jdouble, (Vector3D, Vector3D), arg0, arg1)
end

function distance(obj::Vector3D, arg0::Point)
    return jcall(obj, "distance", jdouble, (Point,), arg0)
end

function distance1(::Type{Vector3D}, arg0::Vector3D, arg1::Vector3D)
    return jcall(Vector3D, "distance1", jdouble, (Vector3D, Vector3D), arg0, arg1)
end

function distance1(obj::Vector3D, arg0::JVector)
    return jcall(obj, "distance1", jdouble, (JVector,), arg0)
end

function distance_inf(::Type{Vector3D}, arg0::Vector3D, arg1::Vector3D)
    return jcall(Vector3D, "distanceInf", jdouble, (Vector3D, Vector3D), arg0, arg1)
end

function distance_inf(obj::Vector3D, arg0::JVector)
    return jcall(obj, "distanceInf", jdouble, (JVector,), arg0)
end

function distance_sq(::Type{Vector3D}, arg0::Vector3D, arg1::Vector3D)
    return jcall(Vector3D, "distanceSq", jdouble, (Vector3D, Vector3D), arg0, arg1)
end

function distance_sq(obj::Vector3D, arg0::JVector)
    return jcall(obj, "distanceSq", jdouble, (JVector,), arg0)
end

function dot_product(::Type{Vector3D}, arg0::Vector3D, arg1::Vector3D)
    return jcall(Vector3D, "dotProduct", jdouble, (Vector3D, Vector3D), arg0, arg1)
end

function dot_product(obj::Vector3D, arg0::JVector)
    return jcall(obj, "dotProduct", jdouble, (JVector,), arg0)
end

function equals(obj::Vector3D, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_alpha(obj::Vector3D)
    return jcall(obj, "getAlpha", jdouble, ())
end

function get_delta(obj::Vector3D)
    return jcall(obj, "getDelta", jdouble, ())
end

function get_norm(obj::Vector3D)
    return jcall(obj, "getNorm", jdouble, ())
end

function get_norm1(obj::Vector3D)
    return jcall(obj, "getNorm1", jdouble, ())
end

function get_norm_inf(obj::Vector3D)
    return jcall(obj, "getNormInf", jdouble, ())
end

function get_norm_sq(obj::Vector3D)
    return jcall(obj, "getNormSq", jdouble, ())
end

function get_space(obj::Vector3D)
    return jcall(obj, "getSpace", Space, ())
end

function get_x(obj::Vector3D)
    return jcall(obj, "getX", jdouble, ())
end

function get_y(obj::Vector3D)
    return jcall(obj, "getY", jdouble, ())
end

function get_z(obj::Vector3D)
    return jcall(obj, "getZ", jdouble, ())
end

function get_zero(obj::Vector3D)
    return jcall(obj, "getZero", Vector3D, ())
end

function hash_code(obj::Vector3D)
    return jcall(obj, "hashCode", jint, ())
end

function is_infinite(obj::Vector3D)
    return jcall(obj, "isInfinite", jboolean, ())
end

function is_na_n(obj::Vector3D)
    return jcall(obj, "isNaN", jboolean, ())
end

function negate(obj::Vector3D)
    return jcall(obj, "negate", Vector3D, ())
end

function normalize(obj::Vector3D)
    return jcall(obj, "normalize", Vector3D, ())
end

function orthogonal(obj::Vector3D)
    return jcall(obj, "orthogonal", Vector3D, ())
end

function scalar_multiply(obj::Vector3D, arg0::jdouble)
    return jcall(obj, "scalarMultiply", Vector3D, (jdouble,), arg0)
end

function subtract(obj::Vector3D, arg0::JVector)
    return jcall(obj, "subtract", Vector3D, (JVector,), arg0)
end

function subtract(obj::Vector3D, arg0::jdouble, arg1::JVector)
    return jcall(obj, "subtract", Vector3D, (jdouble, JVector), arg0, arg1)
end

function to_array(obj::Vector3D)
    return jcall(obj, "toArray", Vector{jdouble}, ())
end

function to_string(obj::Vector3D)
    return jcall(obj, "toString", JString, ())
end

function to_string(obj::Vector3D, arg0::NumberFormat)
    return jcall(obj, "toString", JString, (NumberFormat,), arg0)
end


text:
using Jockey

pwd() == joinpath(@__DIR__, "bin") && cd(@__DIR__) # allow starting app from bin/ dir
push!(Base.modules_warned_for, Base.PkgId(Jockey))
Jockey.main()


text:
using Chain: @chain
using Dates: today, Day

using Stonks: UpdatableSymbol, Symbols, construct_updatable_symbols, build_fx_pair
using Stonks.APIClients: APIClient, APIResource, get_resource, get_type_param
using Stonks.Models:
  AbstractStonksRecord,
  AssetPrice,
  AssetInfo,
  ExchangeRate,
  IncomeStatement,
  BalanceSheet,
  CashflowStatement,
  Earnings
using Stonks.Parsers: AbstractContentParser
using Stonks.Requests: optimistic_request_resolution, prepare_requests, materialize_request

"""
    get_price(
      symbols, [client], [T<:AbstractStonksRecord];
      [interval] = "1d", [from] = missing, [to] = missing, [kwargs...]
    ) -> Union{Vector{T}, Exception}

Retrieves historical time series price data from the configured API client.

### Arguments
- `symbols` can be:
    - `String` with one symbol / ticker
    - `Vector{String}` with multiple symbols
    - `Vector{Tuple{String, Date}}`: tuples of form (symbol, from)
    - `Vector{Tuple{String, Date, Date}}`, tuples of form (symbol, from, to)
- `[client::APIClient]`: can be ommited if one of the correct environmental variable is set (`YAHOOFINANCE_TOKEN` or `ALPHAVANTAGE_TOKEN`)
- `[T<:AbstractStonksRecord]`: data type used for parsing. Change it only if you want to use your custom model. default = `AssetPrice`

### Keywords
- `[interval]`: values = [1d, 1wk, 1mo]. Frequency lower than daily is not supported. default = 1d
- `[from]`: a Date oject indicating lower date limit. default = `missing`
- `[to]`: a Date objject indicating upper date limit. default = `missing`
- `[kwargs...]`: use it to pass keyword arguments if you have url / query parameters that need to be resolved at runtime.

### Examples
```julia-repl
julia> today = Dates.today()
2022-02-18
julia> length(ENV["YAHOOFINANCE_TOKEN"]) # value will be used to create a client
40
julia> get_price("AAPL", from = today - Dates.Day(1))
1-element Vector{AbstractStonksRecord}
 AssetPrice("AAPL", Date("2022-02-17"), 168.88, missing, missing, missing, missing, missing)
 AssetPrice("AAPL", Date("2022-02-18"), 167.3, missing, missing, missing, missing, missing)

julia> from = today - Dates.Day(2)
2022-02-16
julia> prices = get_price(["AAPL", "MSFT"], from = from)
AssetPrice("MSFT", Date("2022-02-16"), 299.5, missing, missing, missing, missing, missing)
AssetPrice("MSFT", Date("2022-02-17"), 290.73, missing, missing, missing, missing, missing
AssetPrice("MSFT", Date("2022-02-18"), 287.93, missing, missing, missing, missing, missing
AssetPrice("AAPL", Date("2022-02-16"), 172.55, missing, missing, missing, missing, missing
AssetPrice("AAPL", Date("2022-02-17"), 168.88, missing, missing, missing, missing, missing
AssetPrice("AAPL", Date("2022-02-18"), 167.3, missing, missing, missing, missing, missing)

julia> prices = get_price([
  ("AAPL", Date("2022-02-17")),
  ("MSFT", Date("2022-02-18"))
  ])
3-element Vector{AbstractStonksRecord}:
 AssetPrice("AAPL", Date("2022-02-17"), 168.88, missing, missing, missing, missing, missing)
 AssetPrice("AAPL", Date("2022-02-18"), 167.3, missing, missing, missing, missing, missing)
 AssetPrice("MSFT", Date("2022-02-18"), 287.93, missing, missing, missing, missing, missing)

julia> prices = get_price([
  ("AAPL", Date("2022-02-15"), Date("2022-02-16")),
  ("MSFT", Date("2022-02-14"), Date("2022-02-15"))
  ])
5-element Vector{AbstractStonksRecord}:
 AssetPrice("MSFT", Date("2022-02-14"), 295.0, missing, missing, missing, missing, missing)
 AssetPrice("MSFT", Date("2022-02-15"), 300.47, missing, missing, missing, missing, missing)
 AssetPrice("AAPL", Date("2022-02-15"), 172.79, missing, missing, missing, missing, missing)
 AssetPrice("AAPL", Date("2022-02-16"), 172.55, missing, missing, missing, missing, missing)
```
"""
function get_price(
  symbols::Symbols,
  client::Union{APIClient,Nothing}=nothing,
  ::Type{T}=AssetPrice;
  interval="1d",
  from=missing,
  to=missing,
  kwargs...,
)::Union{Vector{T},Exception} where {T<:AbstractStonksRecord}
  resource = get_resource(client, T)
  tickers = construct_updatable_symbols(symbols)
  return get_data(resource, tickers; interval=interval, from=from, to=to, kwargs...)
end

"""
    get_info(symbols, [client], [T<:AbstractStonksRecord]) -> Union{Vector{T}, Exception}

Retrieves general information about `symbols`.

### Arguments
- `symbols::Union{String, Vector{String}}`
- `[client]::APIClient` : can be ommited if one of the correct environmental variable is set (`YAHOOFINANCE_TOKEN` or `ALPHAVANTAGE_TOKEN`)
- `[T<:AbstractStonksRecord]` : data type used for parsing. Change it only if you want to use your custom model. default = `AssetInfo`

### Examples
```julia-repl
julia> get_info(client, "AAPL")
1-element Vector{AssetInfo}:
 AssetInfo("AAPL", "USD", "Apple Inc.", "EQUITY", "NMS", "United States",
  "Consumer Electronics", "Technology", "America/New_York", 100000)

julia> length(ENV["YAHOOFINANCE_TOKEN"]) # value will be used to create a client
40
julia> get_info(["AAPL", "MSFT"])
2-element Vector{AssetInfo}:
 AssetInfo("AAPL", "USD", "Apple Inc.", "EQUITY", "NMS", "United States",
  "Consumer Electronics", "Technology", "America/New_York", 100000)
 AssetInfo("MSFT", "USD", "Microsoft Corporation", "EQUITY", "NMS", "United States",
  "Software—Infrastructure", "Technology", "America/New_York", 181000)
```
"""
function get_info(
  symbols::Union{String,Vector{String}},
  client::Union{APIClient,Nothing}=nothing,
  ::Type{T}=AssetInfo;
  kwargs...,
) where {T<:AbstractStonksRecord}
  resource = get_resource(client, T)
  tickers = construct_updatable_symbols(symbols)
  return get_data(resource, tickers, kwargs...)
end

"""
    get_exchange_rate(
      symbols, [client], [T<:AbstractStonksRecord];
      [from] = missing, [to] = missing
    ) -> Union{Vector{T}, Exception}

Retrieves historical exchange rate information

### Arguments
- `symbols` can be:
    - `String` formated as base/quote, each having exactly 3 characters, e.g.: 'EUR/USD', 'USD/CAD'
    - `Vector{String}` with multiple symbols
    - `Vector{Tuple{String, Date}}`: tuples of form (symbol, from)
    - `Vector{Tuple{String, Date, Date}}`, tuples of form (symbol, from, to)
- `[client]::APIClient` can be omitted if one of the correct environmental variable is set (`YAHOOFINANCE_TOKEN` or `ALPHAVANTAGE_TOKEN`)
- `[T<:AbstractStonksRecord]` is the data type used for parsing. Change it only if you want to use your custom model. default = `ExchangeRate`

### Keywords
- `[from]`: a Date object. default = `missing`
- `[to]`: a Date object. default = `missing`

### Examples
```julia-repl
julia> get_exchange_rate("EUR/USD", from=ref_date-Day(1), to=ref_date)
3-element Vector{ExchangeRate}:
 ExchangeRate("EUR", "USD", Date("2022-02-18"), 1.13203)
 ExchangeRate("EUR", "USD", Date("2022-02-17"), 1.13592)
julia> get_exchange_rate(["EUR/USD", "USD/CAD"], from=ref_date-Day(1), to=ref_date)
# 4-element Vector{ExchangeRate}:
 ExchangeRate("EUR", "USD", Date("2022-02-18"), 1.13203)
 ExchangeRate("EUR", "USD", Date("2022-02-17"), 1.13592)
 ExchangeRate("USD", "CAD", Date("2022-02-18"), 1.2748)
 ExchangeRate("USD", "CAD", Date("2022-02-17"), 1.2707)
# Also works with []Tuple{String, Date} or []Tuple{String, Date, Date}
julia>get_exchange_rate([
  ("EUR/USD", Date("2022-02-15"), Date("2022-02-16")),
  ("USD/CAD", Date("2022-02-14"), Date("2022-02-15")),
])
4-element Vector{ExchangeRate}:
...
```
"""
function get_exchange_rate(
  symbols::Symbols,
  client::Union{APIClient,Nothing}=nothing,
  ::Type{T}=ExchangeRate;
  interval="1d",
  from::Union{Date,Missing}=missing,
  to::Union{Date,Missing}=missing,
  kwargs...,
)::Union{Vector{T},Exception} where {T<:AbstractStonksRecord}
  invalid_pairs = @chain symbols begin
    isa(_, String) ? [_] : _
    map(build_fx_pair, _)
    filter(s -> typeof(s) <: Exception, _)
  end
  !isempty(invalid_pairs) && return first(invalid_pairs)
  resource = get_resource(client, T)
  tickers = construct_updatable_symbols(symbols)
  return get_data(resource, tickers; interval=interval, from=from, to=to, kwargs...)
end

"""
    get_income_statement(
      symbols, [client], [T<:AbstractStonksRecord];
      [frequency] = missing, [from] = missing, [to] = missing, [kwargs...]
    ) -> Union{Vector{T}, Exception}

Retrieves main information of income (profit and loss) statement.

### Arguments
- `symbols` can be:
    - `String` with one symbol / ticker
    - `Vector{String}` with multiple symbols
    - `Vector{Tuple{String, Date}}`: tuples of form (symbol, from)
    - `Vector{Tuple{String, Date, Date}}`, tuples of form (symbol, from, to)
- `[client::APIClient]`: can be ommited if one of the correct environmental variable is set (`YAHOOFINANCE_TOKEN` or `ALPHAVANTAGE_TOKEN`)
- `[T<:AbstractStonksRecord]`: data type used for parsing. Change it only if you want to use your custom model. default = `IncomeStatement`

### Keywords
- `[frequency]`: values = [yearly, quarterly]. default = `missing`, which will include both yearly and quarterly frequencies.
- `[from]`: a Date oject indicating lower date limit. default = `missing`
- `[to]`: a Date objject indicating upper date limit. default = `missing`
- `[kwargs...]`: use it to pass keyword arguments if you have url / query parameters that need to be resolved at runtime.

### Examples
```julia
get_income_statement("AAPL")
get_income_statement(["AAPL", "IBM"])
get_income_statement(["AAPL", "IBM"]; frequency="yearly")
get_income_statement([
  ("AAPL", Date("2020-01-01"), Date("2020-12-31")),
  ("MSFT", Date("2020-01-01"), Date("2020-12-31")),
]; frequency="quarterly")
```
"""
function get_income_statement(
  symbols::Symbols,
  client::Union{APIClient,Nothing}=nothing,
  ::Type{T}=IncomeStatement;
  frequency::Union{String,Missing}=missing,
  from::Union{Date,Missing}=missing,
  to::Union{Date,Missing}=missing,
  kwargs...,
)::Union{Vector{T},Exception} where {T<:AbstractStonksRecord}
  resource = get_resource(client, T)
  tickers = construct_updatable_symbols(symbols)
  return get_data(resource, tickers; frequency=frequency, from=from, to=to, kwargs...)
end

"""
    get_balance_sheet(
      symbols, [client], [T<:AbstractStonksRecord];
      [frequency] = missing, [from] = missing, [to] = missing, [kwargs...]
    ) -> Union{Vector{T}, Exception}

Retrieves main information from balance sheet statement.

### Arguments
- `symbols` can be:
    - `String` with one symbol / ticker
    - `Vector{String}` with multiple symbols
    - `Vector{Tuple{String, Date}}`: tuples of form (symbol, from)
    - `Vector{Tuple{String, Date, Date}}`, tuples of form (symbol, from, to)
- `[client::APIClient]`: can be ommited if one of the correct environmental variable is set (`YAHOOFINANCE_TOKEN` or `ALPHAVANTAGE_TOKEN`)
- `[T<:AbstractStonksRecord]`: data type used for parsing. Change it only if you want to use your custom model. default = `BalanceSheet`

### Keywords
- `[frequency]`: values = [yearly, quarterly]. default = `missing`, which will include both yearly and quarterly frequencies.
- `[from]`: a Date oject indicating lower date limit. default = `missing`
- `[to]`: a Date objject indicating upper date limit. default = `missing`
- `[kwargs...]`: use it to pass keyword arguments if you have url / query parameters that need to be resolved at runtime.

### Examples
```julia
get_balance_sheet("AAPL")
get_balance_sheet(["AAPL", "IBM"])
get_balance_sheet(["AAPL", "IBM"]; frequency="yearly")
get_balance_sheet([
  ("AAPL", Date("2020-01-01"), Date("2020-12-31")),
  ("MSFT", Date("2020-01-01"), Date("2020-12-31")),
]; frequency="quarterly")
```
"""
function get_balance_sheet(
  symbols::Symbols,
  client::Union{APIClient,Nothing}=nothing,
  ::Type{T}=BalanceSheet;
  frequency::Union{String,Missing}=missing,
  from::Union{Date,Missing}=missing,
  to::Union{Date,Missing}=missing,
  kwargs...,
)::Union{Vector{T},Exception} where {T<:AbstractStonksRecord}
  resource = get_resource(client, T)
  tickers = construct_updatable_symbols(symbols)
  return get_data(resource, tickers; frequency=frequency, from=from, to=to, kwargs...)
end

"""
    get_cashflow_statement(
      symbols, [client], [T<:AbstractStonksRecord];
      [frequency] = missing, [from] = missing, [to] = missing, [kwargs...]
    ) -> Union{Vector{T}, Exception}

Retrieves main information from cash flow statement.

### Arguments
- `symbols` can be:
    - `String` with one symbol / ticker
    - `Vector{String}` with multiple symbols
    - `Vector{Tuple{String, Date}}`: tuples of form (symbol, from)
    - `Vector{Tuple{String, Date, Date}}`, tuples of form (symbol, from, to)
- `[client::APIClient]`: can be ommited if one of the correct environmental variable is set (`YAHOOFINANCE_TOKEN` or `ALPHAVANTAGE_TOKEN`)
- `[T<:AbstractStonksRecord]`: data type used for parsing. Change it only if you want to use your custom model. default = `CashflowStatement`

### Keywords
- `[frequency]`: values = [yearly, quarterly]. default = `missing`, which will include both yearly and quarterly frequencies.
- `[from]`: a Date oject indicating lower date limit. default = `missing`
- `[to]`: a Date objject indicating upper date limit. default = `missing`
- `[kwargs...]`: use it to pass keyword arguments if you have url / query parameters that need to be resolved at runtime.

### Examples
```julia
get_cashflow_statement("AAPL")
get_cashflow_statement(["AAPL", "IBM"])
get_cashflow_statement(["AAPL", "IBM"]; frequency="yearly")
get_cashflow_statement([
  ("AAPL", Date("2020-01-01"), Date("2020-12-31")),
  ("MSFT", Date("2020-01-01"), Date("2020-12-31")),
]; frequency="quarterly")
```
"""
function get_cashflow_statement(
  symbols::Symbols,
  client::Union{APIClient,Nothing}=nothing,
  ::Type{T}=CashflowStatement;
  frequency::Union{String,Missing}=missing,
  from::Union{Date,Missing}=missing,
  to::Union{Date,Missing}=missing,
  kwargs...,
)::Union{Vector{T},Exception} where {T<:AbstractStonksRecord}
  resource = get_resource(client, T)
  tickers = construct_updatable_symbols(symbols)
  return get_data(resource, tickers; frequency=frequency, from=from, to=to, kwargs...)
end

"""
    get_earnings(
      symbols, [client], [T<:AbstractStonksRecord];
      [frequency] = missing, [from] = missing, [to] = missing, [kwargs...]
    ) -> Union{Vector{T}, Exception}

Retrieves historical informaton about earnings per share announcements.

### Arguments
- `symbols` can be:
    - `String` with one symbol / ticker
    - `Vector{String}` with multiple symbols
    - `Vector{Tuple{String, Date}}`: tuples of form (symbol, from)
    - `Vector{Tuple{String, Date, Date}}`, tuples of form (symbol, from, to)
- `[client::APIClient]`: can be ommited if one of the correct environmental variable is set (`YAHOOFINANCE_TOKEN` or `ALPHAVANTAGE_TOKEN`)
- `[T<:AbstractStonksRecord]`: data type used for parsing. Change it only if you want to use your custom model. default = `Earnings`

### Keywords
- `[frequency]`: values = [yearly, quarterly]. default = `missing`, which will include both yearly and quarterly frequencies.
- `[from]`: a Date oject indicating lower date limit. default = `missing`
- `[to]`: a Date objject indicating upper date limit. default = `missing`
- `[kwargs...]`: use it to pass keyword arguments if you have url / query parameters that need to be resolved at runtime.

### Examples
```julia
get_earnings("AAPL")
get_earnings(["AAPL", "IBM"])
get_earnings(["AAPL", "IBM"]; frequency="yearly")
get_earnings([
  ("AAPL", Date("2020-01-01"), Date("2020-12-31")),
  ("MSFT", Date("2020-01-01"), Date("2020-12-31")),
]; frequency="quarterly")
```
"""
function get_earnings(
  symbols::Symbols,
  client::Union{APIClient,Nothing}=nothing,
  ::Type{T}=Earnings;
  frequency::Union{String,Missing}=missing,
  from::Union{Date,Missing}=missing,
  to::Union{Date,Missing}=missing,
  kwargs...,
)::Union{Vector{T},Exception} where {T<:AbstractStonksRecord}
  resource = get_resource(client, T)
  tickers = construct_updatable_symbols(symbols)
  return get_data(resource, tickers; frequency=frequency, from=from, to=to, kwargs...)
end

"""
    get_data(resource, symbols; kwargs...) -> Union{Vector{<:AbstractStonksRecord}, Exception}

Generic function to get data of type resource{T}. Functions such as get_price, get_info, call this.

### Arguments
- `resource::APIResource`: instance of an `APIResource`.
- `symbols::Symbol` can be:
    - `String` with one symbol / ticker
    - `Vector{String}` with multiple symbols
    - `Vector{Tuple{String, Date}}`: tuples of form (symbol, from)
    - `Vector{Tuple{String, Date, Date}}`, tuples of form (symbol, from, to)
    - `Vector{UpdatableSymbol}`: type Stonks.UpdatableSymbol(symbol; [from], [to])

### Keywords
- `[interval]`: values = 1d, 1wk, 1mo. Frequency lower than daily is not supported. default = 1d
- `[from]`: a Date oject indicating lower date limit. default = `missing`
- `[to]`: a Date objject indicating upper date limit. default = `missing`
- `[kwargs...]`: use it to pass keyword arguments if you have url / query parameters that need to be resolved at runtime.
"""
function get_data(
  resource::Union{APIResource,Exception},
  symbols::Union{Vector{UpdatableSymbol},Symbols,Missing}=missing;
  kwargs...,
)#::Union{Vector{AbstractStonksRecord},Exception}
  tickers = (
    if ismissing(symbols)
      [UpdatableSymbol("FOO")]
    elseif isa(symbols, Vector{UpdatableSymbol})
      symbols
    else
      construct_updatable_symbols(symbols)
    end
  )
  typeof(resource) <: Exception && return resource
  requests = prepare_requests(tickers, resource; kwargs...)
  typeof(requests) <: Exception && return requests
  responses = Channel(length(requests))
  @sync begin
    for rp in requests
      Threads.@spawn begin
        symbol = length(rp.tickers) == 1 ? first(rp.tickers).ticker : missing
        frequency = get(kwargs, :frequency, missing)
        res = materialize_request(
          rp, resource.parser; from=rp.from, to=rp.to, symbol=symbol, frequency=frequency
        )
        err = typeof(res) <: Exception ? res : nothing
        value = typeof(res) <: Exception ? nothing : res
        push!(responses, (params=rp, err=err, value=value))
      end
    end
  end
  close(responses)
  return optimistic_request_resolution(get_type_param(resource), responses)
end

# function get_data(
#   client::APIClient,
#   tickers::Vector{UpdatableSymbol},
# )::Dict{String,Union{Vector{AbstractStonksRecord},Exception}}
#   responses = Channel(length(client.resources))
#   @sync begin
#     for (k, resource) in client.resources
#       Threads.@spawn begin
#         resp = get_data(resource, tickers)
#         push!(responses, (resource = k, value = resp))
#       end
#     end
#   end
#   close(responses)
#   collect_results(responses)
# end

# function collect_results(c::Channel)
#   result = Dict()
#   for item in c
#     has_value = isa(item.value, Vector{AbstractStonksRecord})
#     if haskey(result, item.resource) && has_value
#       append!(result[item.resource], item.value)
#     elseif has_value
#       result[item.resource] = item.value
#     else
#       result[item.resource] = item.value
#     end
#   end
#   return result
# end


text:
const EVAL_FUNCTION_NAME = "evaluate"

function evaluate(f::FlatForest, features::Vector{Float64})
    result = 0.0
    for tree in f.trees
        result += evaluate(tree, features)
    end
    return result
end

type CompiledForest <: Evaluator
    handle::Ptr{Void}
    func::Ptr{Void}

    function CompiledForest(f::Forest)
        function load_shared_object(path::String)
            h = dlopen(path, RTLD_LAZY | RTLD_LOCAL)
            @assert h != C_NULL
            p = dlsym(h, EVAL_FUNCTION_NAME)
            @assert p != C_NULL
            result = new(h, p)
            finalizer(result, r -> dlclose(r.handle))
            return result
        end

        function write_code_to_file(f::Forest)
            tempdir = mktempdir()
            cpp_file = joinpath(tempdir, "tree.cpp")
            stream = open(cpp_file, "w")
            buf = IOBuffer()
            write(stream, "extern \"C\" {\n")
            code_gen(f, stream)
            write(stream, "}\n")
            flush(stream)
            return cpp_file
        end

        function compile_to_shared_object(cpp_file::String)
            tempdir = mktempdir()
            object_file = joinpath(tempdir, "tree.o")
            shared_object_file = joinpath(tempdir, "tree.so")
            run(`clang $cpp_file -c -O3 -o $object_file`)
            run(`clang -shared $object_file -dynamiclib -O3 -o $shared_object_file`)
            return shared_object_file
        end

        return f |>
        write_code_to_file |>
        compile_to_shared_object |>
        load_shared_object
    end
end

function evaluate(c::CompiledForest, features::Vector{Float64})
    return ccall(c.func, Float64, (Ptr{Float64},), features)
end

function code_gen(f::Forest, io::IO)
    for i = 1:length(f.trees)
        write(io, "double evaluate_$i(const double* f) {\n")
        code_gen(f.trees[i], io)
        write(io, "}\n")
    end

    write(io, "double $EVAL_FUNCTION_NAME(const double* f) {\n")
    write(io, "double result = 0.0;\n")
    for i = 1:length(f.trees)
        write(io, "result += evaluate_$i(f);\n")
    end
    write(io, "return result;\n")
    write(io, "}\n")
end

function code_gen(t::Leaf, io::IO)
    value = t.value
    write(io, "return $value;\n")
end

function code_gen(i::Inner, io::IO)
    split_value = i.splitValue
    # Switch to C-style array indexing
    feature = i.feature - 1
    write(io, "if (f[$feature] < $split_value) {\n")
    code_gen(i.left, io)
    write(io, "} else {\n")
    code_gen(i.right, io)
    write(io, "}\n")
end


text:
export randomM

"""
function randomM(m,n;p=2,N=1)

Generate samples for M from uniform distribution with non-repeated columns.

Input:
    m - number of blocks
    n - number of columns
    p - number of categories in each block (default=2)
    N - number of samples (default=1)

Output:
    M - BitArray of size (m*(p-1),n,N), or (m*(p-1),n) if N==1
"""
function randomM(m::Int64,n::Int64;p::Int64=2,N::Int64=1)
    if log(n) > m*log(p)
        error("unique columns in M not possible if n > p^m")
    end

    # store all feasible blocks of M in a matrix B of size ((p-1)*p^n, n)
    B = getB(n,p)

    M = falses(m*(p-1),n,N)
    for idx = 1:N
        # rejection sampling: check uniqueness of columns
        firstFlag = true # (this should be a do-while loop)
        while firstFlag || size(unique(M[:,:,idx],2),2) < n
            firstFlag = false
            # randomly choose m blocks from B
            blockIdB = rand(1:p^n, m)
            for blockIdM = 1:m
                subsB = (blockIdB[blockIdM]-1)*(p-1)+1:blockIdB[blockIdM]*(p-1)
                subsM = (blockIdM-1)*(p-1)+1:blockIdM*(p-1)
                M[subsM,:,idx] = B[subsB,:]
            end
        end
    end

    # output
    if N==1
        M = squeeze(M,3)
    end
    return M
end


text:
module DependencyTree
  using AbstractTrees

  export TreeNode

  mutable struct TreeNode
    token::String
    children::Vector{TreeNode}
    parent::TreeNode
    parent_label::String

    TreeNode(token::String) = new(token, Vector{TreeNode}())
    TreeNode(token::String, parent::TreeNode, parent_label::String) = (node = new(token, Vector{TreeNode}(), parent, parent_label); add_child(parent, node); node)
  end

  function AbstractTrees.children(node::TreeNode)
    node.children
  end

  function add_child(parent::TreeNode, child::TreeNode)
    push!(parent.children, child)
  end

  function find_node(root::TreeNode, token::String)
    if root.token == token
      return root
    end

    for node in root.children
      result = find_node(node, token)
      if result === nothing
        continue
      end

      return result
    end

    nothing
  end
end


text:
#=
Author:= Shreeda Bhat
License:= MIT
UTF-8 =#

a = 1-2+3+4
println(a)
println("-------------------------------------")
b=1/2+3/4
println(b)
println("--------------------------------------")
c = (1/2)+(3/4)
println(c)
println("---------------------------------------")
d = 1+(2/3)-4
println(d)
println("---------------------------------------")
println("*******************************")
e = (1+2/3^2)
println(e)
println("#######################################")


text:
using HydroRefStations

using Test
using Dates: Day

@testset "HRS.jl" begin
    hrs = HRS()
    @testset "HRS()" begin
        nrows, ncols = size(hrs.sites)
        @test nrows > 0 # At least one site.
        @test ncols > 0 # At least one column.
        @test length(hrs.header) > 0 # At least one header line
        @test hrs.header[1][1] != '#'
        @test hrs.header[1][1] != ','
        @test hrs.header[end][1] != '#'
        @test hrs.header[end][end] != '\"'

        @test length(hrs.data_types) > 0 # At least one data type.

        data_types = hrs.data_types["year"]
        @test "annual data" in data_types
        @test "11 year moving average" in data_types
        @test "seasonal data" ∉ data_types

        data_types = hrs.data_types["month"]
        @test "monthly data" in data_types
        @test "monthly boxplot" in data_types
        @test "daily data" ∉ data_types
        @test "seasonal data" ∉ data_types
    end

    @testset "get_data()" begin
        awrc_ids = ["410730", "114001A"]
        data_types = ["monthly data", "seasonal data"]
        for awrc_id in awrc_ids
            for data_type in data_types
                data, header = get_data(hrs, awrc_id, data_type)

                local nrows, ncols = size(data)
                @test nrows > 0 # At least one data point.
                @test ncols > 0 # At least one column.

                @test length(header) > 0 # At least one header line
                @test header[end][1] != '#'
                @test header[end][end] != '\"'
                if data_type == "monthly data"
                    # Ignore the first and last 12 months in a case of incomplete dataset.
                    @test maximum(diff(data[!,"Start Date"][13:end-12])) <= Day(31)
                    @test minimum(diff(data[!,"Start Date"][13:end-12])) >= Day(28)
                end
                if data_type == "seasonal data"
                    @test maximum(diff(data[!,"Start Date"][5:end-4])) <= Day(92)
                    @test minimum(diff(data[!,"Start Date"][5:end-4])) >= Day(89)
                end
            end
            @test_throws ArgumentError get_data(hrs, "invalid ID", "daily data")
            @test_throws ArgumentError get_data(hrs, awrc_id, "invalid data type")
        end

        data, header = get_data(hrs, awrc_ids[1], "SEASONAL DATA")
    end
    @testset "show()" begin
        show_str = repr(hrs)
        @test occursin("Hydrologic Reference Stations", show_str)
        @test occursin("AWRC", show_str)
        @test occursin("Catchment", show_str)
    end
    @testset "close()" begin
        close!(hrs)
        @test isempty(hrs.sites)
        @test isempty(hrs.header)
        @test isempty(hrs.data_types)
    end

end


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
AbstractMultivector.jl defines the AbstractMultivector type and interface
"""

# --- Exports

# Types
export AbstractMultivector

# Properties
import Base.getindex
import LinearAlgebra.norm
export dim, grades, blades, norm

# Functions
import Base.:(-), Base.reverse
export dual

import Base.:(==), Base.isapprox
import Base.iszero, Base.isone

import Base.zero, Base.one
import Base.convert

# --- Types

"""
    AbstractMultivector{<:AbstractFloat}

Supertype for all multivector types.

Interface
=========

Properties
----------

    dim(M::AbstractMultivector)::Int

    grades(M::AbstractMultivector)::Vector{Int}

    blades(M::AbstractMultivector)::Vector{<:AbstractBlade}

    norm(M::AbstractMultivector)::AbstractFloat

    getindex(M::AbstractMultivector, k::Int)::Vector{<:AbstractBlade}

Unary Operations
----------------

    -(M::AbstractMultivector)::AbstractMultivector

    reverse(M::AbstractMultivector)::AbstractMultivector

    dual(M::AbstractMultivector)::AbstractMultivector

Binary Operations
-----------------

    +(M::AbstractMultivector, N::AbstractMultivector)::AbstractMultivector

    -(M::AbstractMultivector, N::AbstractMultivector)::AbstractMultivector

    *(M::AbstractMultivector, N::AbstractMultivector)::AbstractMultivector

    /(M::AbstractMultivector, N::AbstractMultivector)::AbstractMultivector

    wedge(M::AbstractMultivector, N::AbstractMultivector)::AbstractMultivector
    ∧(M::AbstractMultivector, N::AbstractMultivector)::AbstractMultivector

    contract_left(M::AbstractMultivector,
                  N::AbstractMultivector)::AbstractMultivector
    <(M::AbstractMultivector, N::AbstractMultivector)::AbstractMultivector

    contract_right(M::AbstractMultivector,
                   N::AbstractMultivector)::AbstractMultivector
    >(M::AbstractMultivector, N::AbstractMultivector)::AbstractMultivector

    dot(M::AbstractMultivector, N::AbstractMultivector;
        left=true)::AbstractMultivector
    ⋅(M::AbstractMultivector, N::AbstractMultivector;
      left=true)::AbstractMultivector

    project(M::AbstractMultivector, B::AbstractBlade)::AbstractMultivector

Comparison Functions
--------------------

    ==(M::AbstractMultivector, N::AbstractMultivector)::Bool
    isequal(M::AbstractMultivector, N::AbstractMultivector)::Bool

    ≈(M::AbstractMultivector, N::AbstractMultivector)::Bool
    isapprox(M::AbstractMultivector, N::AbstractMultivector)::Bool

Utility Functions
-----------------

    convert(::Type{T}, M::AbstractMultivector) where {T<:AbstractFloat}

Implementation
==============

* The return value of all methods should preserve the precision of its
  AbstractMultivector arguments (when possible).

"""
abstract type AbstractMultivector{T<:AbstractFloat} end

# --- Public functions/methods

"""
    dim(M::AbstractMultivector)::Int

Return dimension of the real space that `M` is embedded within.
"""
function dim end

"""
    grades(M::Multivector; collect=true)::Vector{Int}

Return the grades of the nonzero `k`-vector components of `M`. When `collect`
is `false`, an iterator over the grades is returned.
"""
function grades end

"""
    blades(M::AbstractMultivector)::Vector{<:AbstractBlade}

Return the blades that `M` is composed of.
"""
function blades end

"""
    norm(M::AbstractMultivector)::AbstractFloat

Compute the norm of `M`.
"""
function norm end

"""
    getindex(M::AbstractMultivector, k::Int)::Vector{<:AbstractBlade}

Return the `k`-vector component of `M`.

Notes
=====

* When `M` is a blade, `M[k]` is a Vector containing `M` if the grade of `M`
  is equal to `k`; otherwise, `M[k]` is an empty vector.
"""
function getindex end

"""
    reverse(M::AbstractMultivector)::AbstractMultivector

Compute the reverse of `M`.
"""
function reverse end

"""
    dual(M::AbstractMultivector)::AbstractMultivector

Compute the dual of `M` (relative to the pseudoscalar of the geometric algebra
that `M` is an element of).
"""
function dual end

# ------ Comparison methods

# B::AbstractMultivector, C::AbstractMultivector
isapprox(B::AbstractMultivector, C::AbstractMultivector) = false

# B::AbstractMultivector, x::Real
# x::Real, B::AbstractMultivector
isapprox(M::AbstractMultivector, x::Real) = false
isapprox(x::Real, M::AbstractMultivector) = false

# B::AbstractMultivector, v::Vector
# v::Vector, B::AbstractMultivector
isapprox(M::AbstractMultivector, v::Vector) = false
isapprox(v::Vector, M::AbstractMultivector) = false

iszero(M::AbstractMultivector) = (M === zero(M))

isone(M::AbstractMultivector) = (M === one(M))

# ------ Utility methods

zero(M::AbstractMultivector) = Zero{typeof(norm(M))}()
zero(::Type{<:AbstractMultivector}) = Zero{Float64}()
zero(::Type{<:AbstractMultivector{T}}) where {T<:AbstractFloat} = Zero{T}()

one(M::AbstractMultivector) = One{typeof(norm(M))}()
one(::Type{<:AbstractMultivector}) = One{Float64}()
one(::Type{<:AbstractMultivector{T}}) where {T<:AbstractFloat} = One{T}()

"""
    convert(::Type{T}, M::AbstractMultivector) where {T<:AbstractFloat}

Convert AbstractMultivector to have the floating-point precision of type `T`.
"""
function convert end

# --- Non-exported utility methods

"""
    assert_dim_equal(M::AbstractMultivector, N::AbstractMultivector)
    assert_dim_equal(M::AbstractMultivector, v::Vector)
    assert_dim_equal(v::Vector, M::AbstractMultivector)

Assert that the dimensions of `M` and `N` are the same.
"""
function assert_dim_equal(M::AbstractMultivector, N::AbstractMultivector)
    if dim(M) != dim(N)
        throw(DimensionMismatch("`dim(M)` not equal to `dim(N)`"))
    end
end

function assert_dim_equal(M::AbstractMultivector, v::Vector)
    if dim(M) != length(v)
        throw(DimensionMismatch("`dim(M)` not equal to `length(v)`"))
    end
end
assert_dim_equal(v::Vector, M::AbstractMultivector) = assert_dim_equal(M, v)