text:
export null, isnull, all_ones

null(typ::LLVMType) = Value(API.LLVMConstNull(typ))

all_ones(typ::LLVMType) = Value(API.LLVMConstAllOnes(typ))

isnull(val::Value) = convert(Core.Bool, API.LLVMIsNull(val))

abstract type Constant <: User end

unsafe_destroy!(constant::Constant) = API.LLVMDestroyConstant(constant)

# forward declarations
@checked mutable struct Module
    ref::API.LLVMModuleRef
end
abstract type Instruction <: User end

## data

export ConstantData, PointerNull, UndefValue, ConstantInt, ConstantFP

abstract type ConstantData <: Constant end

@checked struct PointerNull <: ConstantData
    ref::API.LLVMValueRef
end
register(PointerNull, API.LLVMConstantPointerNullValueKind)

PointerNull(typ::PointerType) = PointerNull(API.LLVMConstPointerNull(typ))

@checked struct UndefValue <: ConstantData
    ref::API.LLVMValueRef
end
register(UndefValue, API.LLVMUndefValueValueKind)

UndefValue(typ::LLVMType) = UndefValue(API.LLVMGetUndef(typ))

@checked struct ConstantInt <: ConstantData
    ref::API.LLVMValueRef
end
register(ConstantInt, API.LLVMConstantIntValueKind)

# NOTE: fixed set for dispatch, also because we can't rely on sizeof(T)==width(T)
const WideInteger = Union{Int64, UInt64}
ConstantInt(typ::IntegerType, val::WideInteger, signed=false) =
    ConstantInt(API.LLVMConstInt(typ, reinterpret(Culonglong, val),
                convert(Bool, signed)))
const SmallInteger = Union{Core.Bool, Int8, Int16, Int32, UInt8, UInt16, UInt32}
ConstantInt(typ::IntegerType, val::SmallInteger, signed=false) =
    ConstantInt(typ, convert(Int64, val), signed)

function ConstantInt(typ::IntegerType, val::Integer, signed=false)
    valbits = ceil(Int, log2(abs(val))) + 1 # FIXME: doesn't work for val=0
    numwords = ceil(Int, valbits / 64)
    words = Vector{Culonglong}(undef, numwords)
    for i in 1:numwords
        words[i] = (val >> 64(i-1)) % Culonglong
    end
    return ConstantInt(API.LLVMConstIntOfArbitraryPrecision(typ, numwords, words))
end

# NOTE: fixed set where sizeof(T) does match the numerical width
const SizeableInteger = Union{Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32, UInt64, UInt128}
function ConstantInt(val::T; ctx::Context) where T<:SizeableInteger
    typ = IntType(sizeof(T)*8; ctx)
    return ConstantInt(typ, val, T<:Signed)
end

# Booleans are encoded with a single bit, so we can't use sizeof
ConstantInt(val::Core.Bool; ctx::Context) = ConstantInt(Int1Type(ctx), val ? 1 : 0)

Base.convert(::Type{T}, val::ConstantInt) where {T<:Unsigned} =
    convert(T, API.LLVMConstIntGetZExtValue(val))

Base.convert(::Type{T}, val::ConstantInt) where {T<:Signed} =
    convert(T, API.LLVMConstIntGetSExtValue(val))

# Booleans aren't Signed or Unsigned
Base.convert(::Type{Core.Bool}, val::ConstantInt) = convert(Int, val) != 0

@checked struct ConstantFP <: ConstantData
    ref::API.LLVMValueRef
end
register(ConstantFP, API.LLVMConstantFPValueKind)

ConstantFP(typ::FloatingPointType, val::Real) =
    ConstantFP(API.LLVMConstReal(typ, Cdouble(val)))

ConstantFP(val::Float16; ctx::Context) =
    ConstantFP(HalfType(ctx), val)
ConstantFP(val::Float32; ctx::Context) =
    ConstantFP(FloatType(ctx), val)
ConstantFP(val::Float64; ctx::Context) =
    ConstantFP(DoubleType(ctx), val)

Base.convert(::Type{T}, val::ConstantFP) where {T<:AbstractFloat} =
    convert(T, API.LLVMConstRealGetDouble(val, Ref{API.LLVMBool}()))

# sequential

export ConstantDataSequential, ConstantDataArray, ConstantDataVector

abstract type ConstantDataSequential <: Constant end

@checked struct ConstantDataArray <: ConstantDataSequential
    ref::API.LLVMValueRef
end
register(ConstantDataArray, API.LLVMConstantDataArrayValueKind)

@checked struct ConstantDataVector <: ConstantDataSequential
    ref::API.LLVMValueRef
end
register(ConstantDataVector, API.LLVMConstantDataVectorValueKind)

# aggregate zero

export ConstantAggregateZero

@checked struct ConstantAggregateZero <: ConstantData
    ref::API.LLVMValueRef
end
register(ConstantAggregateZero, API.LLVMConstantAggregateZeroValueKind)

# array interface
# FIXME: can we reuse the ::ConstantArray functionality with ConstantAggregateZero values?
#        probably works fine if we just get rid of the refcheck
Base.eltype(caz::ConstantAggregateZero) = eltype(llvmtype(caz))
Base.size(caz::ConstantAggregateZero) = (0,)
Base.length(caz::ConstantAggregateZero) = 0
Base.axes(caz::ConstantAggregateZero) = (Base.OneTo(0),)
Base.collect(caz::ConstantAggregateZero) = Value[]

## regular aggregate

export ConstantAggregate

abstract type ConstantAggregate <: Constant end

# arrays

export ConstantArray

@checked struct ConstantArray <: ConstantAggregate
    ref::API.LLVMValueRef
end
register(ConstantArray, API.LLVMConstantArrayValueKind)
register(ConstantArray, API.LLVMConstantDataArrayValueKind)

ConstantArrayOrAggregateZero(value) = Value(value)::Union{ConstantArray,ConstantAggregateZero}

# generic constructor taking an array of constants
function ConstantArray(typ::LLVMType, data::AbstractArray{T,N}=T[]) where {T<:Constant,N}
    @assert all(x->x==typ, llvmtype.(data))

    if N == 1
        return ConstantArrayOrAggregateZero(API.LLVMConstArray(typ, Array(data), length(data)))
    end

    ca_vec = map(x->ConstantArray(typ, x), eachslice(data, dims=1))
    ca_typ = llvmtype(first(ca_vec))

    return ConstantArray(API.LLVMConstArray(ca_typ, ca_vec, length(ca_vec)))
end

# shorthands with arrays of plain Julia data
# FIXME: duplicates the ConstantInt/ConstantFP conversion rules
# XXX: X[X(...)] instead of X.(...) because of empty-container inference
ConstantArray(data::AbstractArray{T,N}; ctx::Context) where {T<:Integer,N} =
    ConstantArray(IntType(sizeof(T)*8; ctx), ConstantInt[ConstantInt(x; ctx) for x in data])
ConstantArray(data::AbstractArray{Core.Bool,N}; ctx::Context) where {N} =
    ConstantArray(Int1Type(ctx), ConstantInt[ConstantInt(x; ctx) for x in data])
ConstantArray(data::AbstractArray{Float16,N}; ctx::Context) where {N} =
    ConstantArray(HalfType(ctx), ConstantFP[ConstantFP(x; ctx) for x in data])
ConstantArray(data::AbstractArray{Float32,N}; ctx::Context) where {N} =
    ConstantArray(FloatType(ctx), ConstantFP[ConstantFP(x; ctx) for x in data])
ConstantArray(data::AbstractArray{Float64,N}; ctx::Context) where {N} =
    ConstantArray(DoubleType(ctx), ConstantFP[ConstantFP(x; ctx) for x in data])

# convert back to known array types
function Base.collect(ca::ConstantArray)
    constants = Array{Value}(undef, size(ca))
    for I in CartesianIndices(size(ca))
        @inbounds constants[I] = ca[Tuple(I)...]
    end
    return constants
end

# array interface
Base.eltype(ca::ConstantArray) = eltype(llvmtype(ca))
function Base.size(ca::ConstantArray)
    dims = Int[]
    typ = llvmtype(ca)
    while typ isa ArrayType
        push!(dims, length(typ))
        typ = eltype(typ)
    end
    return Tuple(dims)
end
Base.length(ca::ConstantArray) = prod(size(ca))
Base.axes(ca::ConstantArray) = Base.OneTo.(size(ca))

function Base.getindex(ca::ConstantArray, idx::Integer...)
    # multidimensional arrays are represented by arrays of arrays,
    # which we need to 'peel back' by looking at the operand sets.
    # for the final dimension, we use LLVMGetElementAsConstant
    @boundscheck Base.checkbounds_indices(Base.Bool, axes(ca), idx) ||
        throw(BoundsError(ca, idx))
    I = CartesianIndices(size(ca))[idx...]
    for i in Tuple(I)
        if isempty(operands(ca))
            ca = LLVM.Value(API.LLVMGetElementAsConstant(ca, i-1))
        else
            ca = (Base.@_propagate_inbounds_meta; operands(ca)[i])
        end
    end
    return ca
end

# structs

export ConstantStruct

@checked struct ConstantStruct <: ConstantAggregate
    ref::API.LLVMValueRef
end
register(ConstantStruct, API.LLVMConstantStructValueKind)

ConstantStructOrAggregateZero(value) = Value(value)::Union{ConstantStruct,ConstantAggregateZero}

# anonymous
ConstantStruct(values::Vector{<:Constant}; ctx::Context, packed::Core.Bool=false) =
    ConstantStructOrAggregateZero(API.LLVMConstStructInContext(ctx, values, length(values), convert(Bool, packed)))

# named
ConstantStruct(typ::StructType, values::Vector{<:Constant}) =
    ConstantStructOrAggregateZero(API.LLVMConstNamedStruct(typ, values, length(values)))

# create a ConstantStruct from a Julia object
function ConstantStruct(value::T, name=String(nameof(T)); ctx::Context,
                        anonymous::Core.Bool=false, packed::Core.Bool=false) where {T}
    isbitstype(T) || throw(ArgumentError("Can only create a ConstantStruct from an isbits struct"))
    isprimitivetype(T) && throw(ArgumentError("Cannot create a ConstantStruct from a primitive value"))

    constants = Vector{Constant}()
    for fieldname in fieldnames(T)
        field = getfield(value, fieldname)

        if isa(field, Integer)
            push!(constants, ConstantInt(field; ctx))
        elseif isa(field, AbstractFloat)
            push!(constants, ConstantFP(field; ctx))
        else # TODO: nested structs?
            throw(ArgumentError("only structs with boolean, integer and floating point fields are allowed"))
        end
    end

    if anonymous
        ConstantStruct(constants; ctx, packed)
    elseif haskey(types(ctx), name)
        typ = types(ctx)[name]
        if collect(elements(typ)) != llvmtype.(constants)
            throw(ArgumentError("Cannot create struct $name {$(join(llvmtype.(constants), ", "))} as it is already defined in this context as {$(join(elements(typ), ", "))}."))
        end
        ConstantStruct(typ, constants)
    else
        typ = StructType(name; ctx)
        elements!(typ, llvmtype.(constants))
        ConstantStruct(typ, constants)
    end
end

# vectors

export ConstantVector

@checked struct ConstantVector <: ConstantAggregate
    ref::API.LLVMValueRef
end
register(ConstantVector, API.LLVMConstantVectorValueKind)

## constant expressions

export ConstantExpr,

       const_neg, const_nswneg, const_nuwneg, const_fneg, const_not, const_add,
       const_nswadd, const_nuwadd, const_fadd, const_sub, const_nswsub, const_nuwsub,
       const_fsub, const_mul, const_nswmul, const_nuwmul, const_fmul, const_udiv,
       const_sdiv, const_fdiv, const_urem, const_srem, const_frem, const_and, const_or,
       const_xor, const_icmp, const_fcmp, const_shl, const_lshr, const_ashr, const_gep,
       const_inbounds_gep, const_trunc, const_sext, const_zext, const_fptrunc, const_fpext,
       const_uitofp, const_sitofp, const_fptoui, const_fptosi, const_ptrtoint,
       const_inttoptr, const_bitcast, const_addrspacecast, const_zextorbitcast,
       const_sextorbitcast, const_truncorbitcast, const_pointercast, const_intcast,
       const_fpcast, const_select, const_extractelement, const_insertelement,
       const_shufflevector, const_extractvalue, const_insertvalue

@checked struct ConstantExpr <: Constant
    ref::API.LLVMValueRef
end
register(ConstantExpr, API.LLVMConstantExprValueKind)

opcode(ce::ConstantExpr) = API.LLVMGetConstOpcode(ce)

const_neg(val::Constant) =
    Value(API.LLVMConstNeg(val))

const_nswneg(val::Constant) =
    Value(API.LLVMConstNSWNeg(val))

const_nuwneg(val::Constant) =
    Value(API.LLVMConstNUWNeg(val))

const_fneg(val::Constant) =
    Value(API.LLVMConstFNeg(val))

const_not(val::Constant) =
    Value(API.LLVMConstNot(val))

const_add(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstAdd(lhs, rhs))

const_nswadd(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstNSWAdd(lhs, rhs))

const_nuwadd(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstNUWAdd(lhs, rhs))

const_fadd(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstFAdd(lhs, rhs))

const_sub(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstSub(lhs, rhs))

const_nswsub(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstNSWSub(lhs, rhs))

const_nuwsub(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstNUWSub(lhs, rhs))

const_fsub(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstFSub(lhs, rhs))

const_mul(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstMul(lhs, rhs))

const_nswmul(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstNSWMul(lhs, rhs))

const_nuwmul(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstNUWMul(lhs, rhs))

const_fmul(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstFMul(lhs, rhs))

const_udiv(lhs::Constant, rhs::Constant; exact::Base.Bool=false) =
    Value(exact ? API.LLVMConstExactUDiv(lhs, rhs) : API.LLVMConstUDiv(lhs, rhs))

const_sdiv(lhs::Constant, rhs::Constant; exact::Base.Bool=false) =
    Value(exact ? API.LLVMConstExactSDiv(lhs, rhs) : API.LLVMConstSDiv(lhs, rhs))

const_fdiv(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstFDiv(lhs, rhs))

const_urem(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstURem(lhs, rhs))

const_srem(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstSRem(lhs, rhs))

const_frem(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstFRem(lhs, rhs))

const_and(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstAnd(lhs, rhs))

const_or(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstOr(lhs, rhs))

const_xor(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstXor(lhs, rhs))

const_icmp(Predicate::API.LLVMIntPredicate, lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstICmp(Predicate, lhs, rhs))

const_fcmp(Predicate::API.LLVMRealPredicate, lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstFCmp(Predicate, lhs, rhs))

const_shl(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstShl(lhs, rhs))

const_lshr(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstLShr(lhs, rhs))

const_ashr(lhs::Constant, rhs::Constant) =
    Value(API.LLVMConstAShr(lhs, rhs))

const_gep(val::Constant, Indices::Vector{<:Constant}) =
   Value(API.LLVMConstGEP(val, Indices, length(Indices)))

const_inbounds_gep(val::Constant, Indices::Vector{<:Constant}) =
   Value(API.LLVMConstInBoundsGEP(val, Indices, length(indices)))

const_trunc(val::Constant, ToType::LLVMType) =
    Value(API.LLVMConstTrunc(val, ToType))

const_sext(val::Constant, ToType::LLVMType) =
    Value(API.LLVMConstSExt(val, ToType))

const_zext(val::Constant, ToType::LLVMType) =
    Value(API.LLVMConstZExt(val, ToType))

const_fptrunc(val::Constant, ToType::LLVMType) =
    Value(API.LLVMConstFPTrunc(val, ToType))

const_fpext(val::Constant, ToType::LLVMType) =
    Value(API.LLVMConstFPExt(val, ToType))

const_uitofp(val::Constant, ToType::LLVMType) =
    Value(API.LLVMConstUIToFP(val, ToType))

const_sitofp(val::Constant, ToType::LLVMType) =
    Value(API.LLVMConstSIToFP(val, ToType))

const_fptoui(val::Constant, ToType::LLVMType) =
    Value(API.LLVMConstFPToUI(val, ToType))

const_fptosi(val::Constant, ToType::LLVMType) =
    Value(API.LLVMConstFPToSI(val, ToType))

const_ptrtoint(val::Constant, ToType::LLVMType) =
    Value(API.LLVMConstPtrToInt(val, ToType))

const_inttoptr(val::Constant, ToType::LLVMType) =
    Value(API.LLVMConstIntToPtr(val, ToType))

const_bitcast(val::Constant, ToType::LLVMType) =
    Value(API.LLVMConstBitCast(val, ToType))

const_addrspacecast(val::Constant, ToType::LLVMType) =
    Value(API.LLVMConstAddrSpaceCast(val, ToType))

const_zextorbitcast(val::Constant, ToType::LLVMType) =
    Value(API.LLVMConstZExtOrBitCast(val, ToType))

const_sextorbitcast(val::Constant, ToType::LLVMType) =
    Value(API.LLVMConstSExtOrBitCast(val, ToType))

const_truncorbitcast(val::Constant, ToType::LLVMType) =
    Value(API.LLVMConstTruncOrBitCast(val, ToType))

const_pointercast(val::Constant, ToType::LLVMType) =
    Value(API.LLVMConstPointerCast(val, ToType))

const_intcast(val::Constant, ToType::LLVMType, isSigned::Base.Bool) =
    Value(API.LLVMConstIntCast(val, ToType, isSigned))

const_fpcast(val::Constant, ToType::LLVMType) =
    Value(API.LLVMConstFPCast(val, ToType))

const_select(cond::Constant, if_true::Value, if_false::Value) =
    Value(API.LLVMConstSelect(cond, if_true, if_false))

const_extractelement(vector::Constant, index::Constant) =
    Value(API.LLVMConstExtractElement(vector ,index))

const_insertelement(vector::Constant, element::Value, index::Constant) =
    Value(API.LLVMConstInsertElement(vector ,element, index))

const_shufflevector(vector1::Constant, vector2::Constant, mask::Constant) =
    Value(API.LLVMConstShuffleVector(vector1, vector2, mask))

const_extractvalue(agg::Constant, Idx::Vector{<:Integer}) =
   Value(API.LLVMConstExtractValue(agg, Idx, length(Idx)))

const_insertvalue(agg::Constant, element::Constant, Idx::Vector{<:Integer}) =
   Value(API.LLVMConstInsertValue(agg, element, Idx, length(Idx)))

# TODO: alignof, sizeof, block_address

## inline assembly

export InlineAsm

@checked struct InlineAsm <: Constant
    ref::API.LLVMValueRef
end
register(InlineAsm, API.LLVMInlineAsmValueKind)

InlineAsm(typ::FunctionType, asm::String, constraints::String,
          side_effects::Core.Bool, align_stack::Core.Bool=false) =
    InlineAsm(API.LLVMConstInlineAsm(typ, asm, constraints,
                                     convert(Bool, side_effects),
                                     convert(Bool, align_stack)))

## global values

abstract type GlobalValue <: Constant end

export GlobalValue,
       isdeclaration,
       linkage, linkage!,
       section, section!,
       visibility, visibility!,
       dllstorage, dllstorage!,
       unnamed_addr, unnamed_addr!,
       alignment, alignment!

parent(val::GlobalValue) = Module(API.LLVMGetGlobalParent(val))

isdeclaration(val::GlobalValue) = convert(Core.Bool, API.LLVMIsDeclaration(val))

linkage(val::GlobalValue) = API.LLVMGetLinkage(val)
linkage!(val::GlobalValue, linkage::API.LLVMLinkage) =
    API.LLVMSetLinkage(val, linkage)

function section(val::GlobalValue)
  #=
  The following started to fail on LLVM 4.0:
    Context() do ctx
      LLVM.Module("SomeModule"; ctx) do mod
        st = LLVM.StructType("SomeType"; ctx)
        ft = LLVM.FunctionType(st, [st])
        fn = LLVM.Function(mod, "SomeFunction", ft)
        section(fn) == ""
      end
      end
  =#
  section_ptr = API.LLVMGetSection(val)
  return section_ptr != C_NULL ? unsafe_string(section_ptr) : ""
end
section!(val::GlobalValue, sec::String) = API.LLVMSetSection(val, sec)

visibility(val::GlobalValue) = API.LLVMGetVisibility(val)
visibility!(val::GlobalValue, viz::API.LLVMVisibility) =
    API.LLVMSetVisibility(val, viz)

dllstorage(val::GlobalValue) = API.LLVMGetDLLStorageClass(val)
dllstorage!(val::GlobalValue, storage::API.LLVMDLLStorageClass) =
    API.LLVMSetDLLStorageClass(val, storage)

unnamed_addr(val::GlobalValue) = convert(Core.Bool, API.LLVMHasUnnamedAddr(val))
unnamed_addr!(val::GlobalValue, flag::Core.Bool) =
    API.LLVMSetUnnamedAddr(val, convert(Bool, flag))

const AlignedValue = Union{GlobalValue,Instruction}   # load, store, alloca
alignment(val::AlignedValue) = API.LLVMGetAlignment(val)
alignment!(val::AlignedValue, bytes::Integer) = API.LLVMSetAlignment(val, bytes)

## global variables

abstract type GlobalObject <: GlobalValue end

export GlobalVariable, unsafe_delete!,
       initializer, initializer!,
       isthreadlocal, threadlocal!,
       threadlocalmode, threadlocalmode!,
       isconstant, constant!,
       isextinit, extinit!

@checked struct GlobalVariable <: GlobalObject
    ref::API.LLVMValueRef
end
register(GlobalVariable, API.LLVMGlobalVariableValueKind)

GlobalVariable(mod::Module, typ::LLVMType, name::String) =
    GlobalVariable(API.LLVMAddGlobal(mod, typ, name))

GlobalVariable(mod::Module, typ::LLVMType, name::String, addrspace::Integer) =
    GlobalVariable(API.LLVMAddGlobalInAddressSpace(mod, typ,
                                                   name, addrspace))

unsafe_delete!(::Module, gv::GlobalVariable) = API.LLVMDeleteGlobal(gv)

function initializer(gv::GlobalVariable)
    init = API.LLVMGetInitializer(gv)
    init == C_NULL ? nothing : Value(init)
end
initializer!(gv::GlobalVariable, val::Constant) =
  API.LLVMExtraSetInitializer(gv, val)
initializer!(gv::GlobalVariable, ::Nothing) =
  API.LLVMExtraSetInitializer(gv, C_NULL)

isthreadlocal(gv::GlobalVariable) = convert(Core.Bool, API.LLVMIsThreadLocal(gv))
threadlocal!(gv::GlobalVariable, bool) =
  API.LLVMSetThreadLocal(gv, convert(Bool, bool))

isconstant(gv::GlobalVariable) = convert(Core.Bool, API.LLVMIsGlobalConstant(gv))
constant!(gv::GlobalVariable, bool) =
  API.LLVMSetGlobalConstant(gv, convert(Bool, bool))

threadlocalmode(gv::GlobalVariable) = API.LLVMGetThreadLocalMode(gv)
threadlocalmode!(gv::GlobalVariable, mode) =
  API.LLVMSetThreadLocalMode(gv, mode)

isextinit(gv::GlobalVariable) =
  convert(Core.Bool, API.LLVMIsExternallyInitialized(gv))
extinit!(gv::GlobalVariable, bool) =
  API.LLVMSetExternallyInitialized(gv, convert(Bool, bool))


text:
using DFControl, Test
using Plots
test_bands = qe_read_bands_file(joinpath(@__DIR__, "testassets", "test_job", "nscf.out"))
@test plot(test_bands[1], fermi=3)         != nothing
@test plot(test_bands[3],:relative_cart)  != nothing
@test plot(test_bands[4],:relative_cryst) != nothing
@test plot(test_bands)                    != nothing

t_fermi = qe_read_fermi_from_output(joinpath(@__DIR__, "testassets", "test_job", "nscf.out"))
t_eigval = test_bands[1].eigvals[1]
@test DFControl.apply_fermi_level.(test_bands, 3.2)[1].eigvals[1] == t_eigval - 3.2
@test DFControl.apply_fermi_level.(test_bands, joinpath(@__DIR__, "testassets", "test_job", "nscf.out"))[1].eigvals[1] == t_eigval - t_fermi
DFControl.apply_fermi_level!.(test_bands, 3.2)
@test test_bands[1].eigvals[1] == t_eigval - 3.2
t_eigval = test_bands[1].eigvals[1]
DFControl.apply_fermi_level!.(test_bands, joinpath(@__DIR__, "testassets", "test_job", "nscf.out"))
@test test_bands[1].eigvals[1] == t_eigval - t_fermi

ext: jl


text:
using BinDeps
@BinDeps.setup

librsb = library_dependency("librsb")

version = "1.2.0-rc7"
provides(Sources,
    Dict(URI("https://downloads.sourceforge.net/project/librsb/librsb-$version.tar.gz") => librsb))

provides(BuildProcess,
    Dict(Autotools(libtarget="librsb.la",
                   configure_options=["CC=gcc", "FC=gfortran"]) => librsb))

provides(AptGet,
    Dict("librsb-dev" => librsb))

provides(Pacman,
    Dict("librsb" => librsb))

@BinDeps.install Dict(:librsb => :librsb)


text:
@testset "Decomposition" begin
    d = loadgraph(joinpath(testdir, "testdata", "graph-decomposition.jgz"))

    for g in testgraphs(d)
        corenum = @inferred(core_number(g))
        @test corenum == [3, 3, 3, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 0]

        @test @inferred(k_core(g)) == k_core(g, corenum=corenum) == [1:8;]
        @test @inferred(k_core(g, 2)) == k_core(g, 2, corenum=corenum) == [1:16;]
        @test length(k_core(g, 4)) == 0

        @test @inferred(k_shell(g)) == k_shell(g, corenum=corenum) == [1:8;]
        @test @inferred(k_shell(g, 2)) == k_shell(g, 2, corenum=corenum) == [9:16;]
        @test length(k_shell(g, 4)) == 0

        @test @inferred(k_crust(g)) == k_crust(g, corenum=corenum) == [9:21;]
        @test @inferred(k_crust(g, 2)) == k_crust(g, 2, corenum=corenum) == [9:21;]
        @test @inferred(k_crust(g, 4, corenum=corenum)) == [1:21;]

        @test @inferred(k_corona(g, 1)) == k_corona(g, 1, corenum=corenum) == [17:20;]
        @test @inferred(k_corona(g, 2)) == [10, 12, 13, 14, 15, 16]

        add_edge!(g, 1, 1)
        @test_throws ArgumentError k_core(g)
        @test_throws ArgumentError k_shell(g)
        @test_throws ArgumentError k_crust(g)
        @test_throws ArgumentError k_corona(g, 1)
    end
end


text:
using DataFrames
using MixedModels
using MixedModels: dataset
using MixedModelsExtras
using Statistics
using Test

@testset "LMM" begin
    model = fit(MixedModel, @formula(reaction ~ 1 + (1|subj)), dataset(:sleepstudy))
    @test icc(model, :subj) == icc(model, [:subj]) == icc(model)
    @test icc(model, :subj) ≈ 0.37918288298942798287
end

@testset "Binomial" begin
    cbpp = dataset(:cbpp)
    model = fit(MixedModel, @formula((incid/hsz) ~ 1 + (1|herd)),
                cbpp, Binomial(), wts=float(cbpp.hsz))
    @test icc(model, :herd) == icc(model, [:herd]) == icc(model)
    @test icc(model, :herd) ≈ 0.1668 atol=0.0005
end

@testset "Bernoulli" begin
    contra = dataset(:contra)
    modelbern = fit(MixedModel, @formula(use ~ 1 + (1|urban&dist)),
                    contra, Bernoulli(), fast=true)
    # force treating as a Binomial model
    modelbin = fit(MixedModel, @formula(use ~ 1 + (1|urban&dist)),
                   contra, Binomial(), fast=true, wts=ones(length(contra.dist)))
    # Bernoullis are a special case of binomial, so make sure they give the same answer
    @test icc(modelbern, Symbol("urban & dist")) ≈ icc(modelbin, Symbol("urban & dist"))
end

@testset "Poisson" begin
    grouseticks = DataFrame(dataset(:grouseticks))
    grouseticks.ch = grouseticks.height .- mean(grouseticks.height)
    model = fit(MixedModel,
                @formula(ticks ~ 1 + year + ch + (1|index) + (1|brood) + (1|location)),
                grouseticks, Poisson(), fast=true)
    @test icc(model, :index) ≈ 0.13467352262090606652 atol=0.0005
    @test icc(model, [:index, :brood]) ≈ 0.3878770599741494518 atol=0.0005
    @test icc(model, [:index, :brood, :location]) ≈ 0.53293244949745322003 atol=0.0005
    @test icc(model, [:index, :brood, :location]) == icc(model)
end


text:
function polydivisible(x::Int)::Bool
  accu = 0
  for (i,digit) in x |> digits |> reverse |> enumerate
    accu = muladd( accu, 10, digit )
    if accu % i != 0 return false end
  end
  return true
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
###################################################################
##  macro and functions for sfmodel_spec(), _init(), and _fit()  ##
###################################################################

#? ---- macros for sfmodel_spec() and functions for sfmodel_init() ------------

#=
macro sfdist(arg::Vararg)

    (length(arg)==1)||throw("@sfdist in sfmodel_spec should have one and only one word.")

     if arg[1] ∈ (:Trun, :Truncated, :truncated, :trun, :t)
        ha = :t
     elseif arg[1] ∈ (:Half, :half, :h)
        ha = :h
     elseif arg[1] ∈ (:Expo, :expo, :e)
        ha = :e
     elseif arg[1] ∈ (:Trun_Scale, :Trun_Scaling, :trun_scale, :trun_scaling, :ts, :s)
        ha = :s
     elseif arg[1] ∈ (:Half_Scale, :Half_Scaling, :half_scale, :half_scaling, :hs)
        throw("The half-normal model is naturally a scaling-property model. Please use @sfdist(half) specification.")
     else
        throw("The keyword of @sfdist in sfmodel_spec() is specified incorrectly.")
     end

    return (:dist, [ha])
 end
 =#

"""
 sfdist(arg::Vararg)

An argument in `sfmodel_sepc()`. Specify the distribution assumption of the one-sided stochastic
variable (aka inefficiency term) of the model. Possible choices include `truncated` (or `trun`, `t`), `half` (or `h`),
`exponential` (or `expo`, `e`), and `trun_scale` (or `trun_scaling`, `ts`).

See the help on `sfmodel_spec()` for more information.

# Examples
```julia-repl
sfmodel_spec(sfdist(t), ...)
sfmodel_spec(sfdist(h), ...)
```
"""
function sfdist(arg::Vararg)

    (length(arg)==1) || throw("`sfdist()` in sfmodel_spec should have one and only one word.")

     if arg[1] ∈ (Trun, trun, t)
        ha = :t
     elseif arg[1] ∈ (Half, half, h)
        ha = :h
     elseif arg[1] ∈ (Expo, expo, e)
        ha = :e
     elseif arg[1] ∈ (Trun_Scale, trun_scale, trun_scaling, s)
        ha = :s
     elseif arg[1] ∈ (Half_Scale, half_scale, half_scaling)
        throw("The half-normal model is naturally a scaling-property model. Use `sfdist(half)` specification.")
     else
        throw("The inefficiency distribution in `sfdist()` of sfmodel_spec() is specified incorrectly.")
     end

    return (:dist, [ha])
 end

  # -------------------

#=
macro sftype(arg::Vararg)
  (length(arg)==1) || throw("`@sftype()` in sfmodel_spec should be either production, prod, or cost.")
  if (arg[1]!=:production) && (arg[1]!=:prod) && (arg[1]!=:cost)
       throw("@sftype in sfmodel_spec should be either production, prod, or cost.")
  end
  return (:type, collect(:($(arg))))
end
=#

"""
 sftype(arg::Vararg)

An argument in `sfmodel_sepc()`. Specify whether the model is a `production` (or `prod`) frontier
or a `cost` frontier.

See the help on `sfmodel_spec()` for more information.

# Examples
```julia-repl
sftype(production)
sftype(cost)
```
"""
function sftype(arg::Vararg)
   (length(arg)==1) || throw("`sftype()` in sfmodel_spec should be either production, prod, or cost.")
   if (arg[1]!=production) && (arg[1]!=prod) && (arg[1]!=cost)
        throw("`sftype()` in sfmodel_spec should be either production, prod, or cost.")
   end
   return (:type, [Symbol(arg[1])])
end

  # -------------------
 #=
 macro model(arg::Vararg)
     (arg[1]==:BC1995) ||(arg[1]==:Wang2002)  || throw("@model only supports BC1995, Wang2002")
     return (:type, collect(:($(arg))))
 end
  =#

"""
@depvar(arg::Vararg)

An argument in `sfmodel_sepc()`. Specify the dependent variable using a
column name from a DataFrame.

See the help on `sfmodel_spec()` for more information.

# Examples
```julia-repl
julia> df
100×5 DataFrame
│ Row │ yvar  │ xvar1     │ xvar2     │ zvar      │ _cons   │
│     │ Int64 │ Float64   │ Float64   │ Float64   │ Float64 │
├─────┼───────┼───────────┼───────────┼───────────┼─────────┤
│ 1   │ 1     │ 0.0306449 │ 0.452148  │ 0.808817  │ 1.0     │
│ 2   │ 2     │ 0.460691  │ 0.296092  │ 0.454545  │ 1.0     │
│ 3   │ 3     │ 0.897503  │ 0.376972  │ 0.907454  │ 1.0     │
│ 4   │ 4     │ 0.682894  │ 0.776861  │ 0.161721  │ 1.0     │
⋮
│ 96  │ 96    │ 0.329647  │ 0.0914057 │ 0.825032  │ 1.0     │
│ 97  │ 97    │ 0.0781165 │ 0.338999  │ 0.761652  │ 1.0     │
│ 98  │ 98    │ 0.41394   │ 0.0063118 │ 0.295372  │ 1.0     │
│ 99  │ 99    │ 0.516381  │ 0.285415  │ 1.91995   │ 1.0     │
│ 100 │ 100   │ 0.944     │ 0.702226  │ -0.539848 │ 1.0     │

sfmodel_spec(@depvar(yvar), ...)
```
"""
 macro depvar(arg::Vararg)
     return (:depvar, collect(:($(arg))))
 end

"""
depvar(arg::Vararg)

An argument in `sfmodel_sepc()`. Specify the maxrix name where the matrix
is used as the data of the dependent variable.

See the help on `sfmodel_spec()` for more information.

# Examples
```julia-repl
julia> ymat
100×1 Matrix{Float64}:
  0.0005852467656204083
 -0.9128698116696892
 -1.1087862468093093
 -1.0714018769972091
 -0.606943663382492
  0.37648407866805467
 -0.1281971631844683
  ⋮
  0.8541895741866585
  1.5109216952026845
 -0.3519833126683764
 -1.0378799750720447
 -0.9990384371507885
  0.18858962788775305

sfmodel_spec(depvar(ymat), ...)
```
"""
 function depvar(arg::Vararg)
    return (:depvar, collect(:($(arg))))
end

  # -------------------
"""
@frontier(arg::Vararg)

An argument in `sfmodel_sepc()`. Specify the variables in the `frontier`
function using column names from a DataFrame.

See the help on `sfmodel_spec()` for more information.

# Examples
```julia-repl
julia> df
100×5 DataFrame
│ Row │ yvar  │ xvar1     │ xvar2     │ zvar      │ _cons   │
│     │ Int64 │ Float64   │ Float64   │ Float64   │ Float64 │
├─────┼───────┼───────────┼───────────┼───────────┼─────────┤
│ 1   │ 1     │ 0.0306449 │ 0.452148  │ 0.808817  │ 1.0     │
│ 2   │ 2     │ 0.460691  │ 0.296092  │ 0.454545  │ 1.0     │
│ 3   │ 3     │ 0.897503  │ 0.376972  │ 0.907454  │ 1.0     │
│ 4   │ 4     │ 0.682894  │ 0.776861  │ 0.161721  │ 1.0     │
⋮
│ 96  │ 96    │ 0.329647  │ 0.0914057 │ 0.825032  │ 1.0     │
│ 97  │ 97    │ 0.0781165 │ 0.338999  │ 0.761652  │ 1.0     │
│ 98  │ 98    │ 0.41394   │ 0.0063118 │ 0.295372  │ 1.0     │
│ 99  │ 99    │ 0.516381  │ 0.285415  │ 1.91995   │ 1.0     │
│ 100 │ 100   │ 0.944     │ 0.702226  │ -0.539848 │ 1.0     │

sfmodel_spec(@frontier(xvar1, xvar2, _cons), ...)
```
"""
 macro frontier(arg::Vararg)   # for sfmodel_spec()
     return (:frontier, collect(:($(arg))))
 end

"""
frontier(arg::Vararg)

An argument in `sfmodel_sepc()`. Specify the name of the maxrix
used as the data of the `frontier` function.

See the help on `sfmodel_spec()` for more information.

# Examples
```julia-repl
julia> xmat
100×3 Matrix{Float64}:
 -0.682942    0.440045   1.0
 -0.680902   -1.68124    1.0
  1.29108    -1.5516     1.0
  0.683652   -0.0319451  1.0
 -0.973079    1.11636    1.0
 -0.343229    0.314457   1.0
  0.107583    0.688177   1.0
  ⋮
  0.0943377  -0.781928   1.0
 -0.599142   -1.01591    1.0
 -0.56726    -1.03394    1.0
  1.33522     0.135763   1.0
  1.13235     0.0177493  1.0
 -0.310638   -0.314166   1.0

sfmodel_spec(frontier(xmat), ...)
```
"""
 function frontier(arg::Vararg)    # for sfmodel_init(0.1, 0.1)
     return (:frontier, collect(:($(arg))) )
 end

"""
frontier(arg::Vector)

An argument in `sfmodel_init()`. Specify the initial values for coefficients in
the `frontier` function.

See the help on `sfmodel_init()` for more information.

# Examples
```julia-repl
sfmodel_init(frontier(0.1, 0.2, 0.5), ...)
b0 = ones(3)*0.1
sfmodel_init( frontier(b0), ...)
```
"""
 function frontier(arg::Vector)    # for smodel_init(b0)
     return (:frontier,  arg)
 end

  # -------------------

  """
  @μ(arg::Vararg)

  An argument in `sfmodel_sepc()`. Specify the variables in the `μ`
  function using column names from a DataFrame.

  See the help on `sfmodel_spec()` for more information.

  # Examples
  ```julia-repl
  julia> df
  100×5 DataFrame
  │ Row │ yvar  │ xvar1     │ xvar2     │ zvar      │ _cons   │
  │     │ Int64 │ Float64   │ Float64   │ Float64   │ Float64 │
  ├─────┼───────┼───────────┼───────────┼───────────┼─────────┤
  │ 1   │ 1     │ 0.0306449 │ 0.452148  │ 0.808817  │ 1.0     │
  │ 2   │ 2     │ 0.460691  │ 0.296092  │ 0.454545  │ 1.0     │
  │ 3   │ 3     │ 0.897503  │ 0.376972  │ 0.907454  │ 1.0     │
  │ 4   │ 4     │ 0.682894  │ 0.776861  │ 0.161721  │ 1.0     │
  ⋮
  │ 96  │ 96    │ 0.329647  │ 0.0914057 │ 0.825032  │ 1.0     │
  │ 97  │ 97    │ 0.0781165 │ 0.338999  │ 0.761652  │ 1.0     │
  │ 98  │ 98    │ 0.41394   │ 0.0063118 │ 0.295372  │ 1.0     │
  │ 99  │ 99    │ 0.516381  │ 0.285415  │ 1.91995   │ 1.0     │
  │ 100 │ 100   │ 0.944     │ 0.702226  │ -0.539848 │ 1.0     │

  sfmodel_spec( @μ(zvar, _cons), ...)
  ```
  """
 macro μ(arg::Vararg)
     return (:μ, collect(:($(arg))))
 end

 """
 μ(arg::Vararg)

 An argument in `sfmodel_sepc()`. Specify the name of the maxrix
 used as the data of the `μ` function.

 See the help on `sfmodel_spec()` for more information.

 # Examples
 ```julia-repl
 julia> mumat
 100×2 Matrix{Float64}:
  -0.943133   1.0
  -0.897392   1.0
   0.585447   1.0
  -0.46106    1.0
  -0.54563    1.0
  -0.619428   1.0
   0.0575559  1.0
   ⋮
   0.0844192  1.0
  -1.3339     1.0
   1.29332    1.0
   0.691466   1.0
   0.422962   1.0
   0.374425   1.0

 sfmodel_spec(μ(mumat), ...)
 ```
 """
 function μ(arg::Vararg)
     return (:μ, collect(:($(arg))))
 end

"""
 μ(arg::Vector)

An argument in `sfmodel_init()`. Specify the initial values for coefficients in
the `μ` function.

See the help on `sfmodel_init()` for more information.

# Examples
```julia-repl
sfmodel_init( μ(0.1, 0.2, 0.5), ...)

b0 = ones(3)*0.1
sfmodel_init( μ(b0), ...)
```
"""
 function μ(arg::Vector)
     return (:μ, arg)
 end

  # -------------------

  """
  @mu(arg::Vararg)

  alias of @μ. See help on @μ.
  """
 macro mu(arg::Vararg)
     return (:μ, collect(:($(arg))))
 end

 """
 mu(arg::Vararg)

 alias of μ. See help on μ.
 """
 function mu(arg::Vararg)
     return (:μ, collect(:($(arg))))
 end

 """
 mu(arg::Vector)

 alias of μ. See help on μ.
 """
 function mu(arg::Vector)
     return (:μ, arg)
 end

  # -------------------

  """
  @σᵤ²(arg::Vararg)

  An argument in `sfmodel_sepc()`. Specify the variables in the `σᵤ²`
  function using column names from a DataFrame.

  See the help on `sfmodel_spec()` for more information.

  # Examples
  ```julia-repl
  julia> df
  100×5 DataFrame
  │ Row │ yvar  │ xvar1     │ xvar2     │ zvar      │ _cons   │
  │     │ Int64 │ Float64   │ Float64   │ Float64   │ Float64 │
  ├─────┼───────┼───────────┼───────────┼───────────┼─────────┤
  │ 1   │ 1     │ 0.0306449 │ 0.452148  │ 0.808817  │ 1.0     │
  │ 2   │ 2     │ 0.460691  │ 0.296092  │ 0.454545  │ 1.0     │
  │ 3   │ 3     │ 0.897503  │ 0.376972  │ 0.907454  │ 1.0     │
  │ 4   │ 4     │ 0.682894  │ 0.776861  │ 0.161721  │ 1.0     │
  ⋮
  │ 96  │ 96    │ 0.329647  │ 0.0914057 │ 0.825032  │ 1.0     │
  │ 97  │ 97    │ 0.0781165 │ 0.338999  │ 0.761652  │ 1.0     │
  │ 98  │ 98    │ 0.41394   │ 0.0063118 │ 0.295372  │ 1.0     │
  │ 99  │ 99    │ 0.516381  │ 0.285415  │ 1.91995   │ 1.0     │
  │ 100 │ 100   │ 0.944     │ 0.702226  │ -0.539848 │ 1.0     │

  sfmodel_spec( @σᵤ²(zvar, _cons), ...)
  ```
  """
 macro σᵤ²(arg::Vararg)
     return (:σᵤ², collect(:($(arg))))
 end

 """
 σᵤ²(arg::Vararg)

 An argument in `sfmodel_sepc()`. Specify the name of the maxrix
 used as the data of the `σᵤ²` function.

 See the help on `sfmodel_spec()` for more information.

 # Examples
 ```julia-repl
 julia> zmat
 100×2 Matrix{Float64}:
  -0.943133   1.0
  -0.897392   1.0
   0.585447   1.0
  -0.46106    1.0
  -0.54563    1.0
  -0.619428   1.0
   0.0575559  1.0
   ⋮
   0.0844192  1.0
  -1.3339     1.0
   1.29332    1.0
   0.691466   1.0
   0.422962   1.0
   0.374425   1.0

 sfmodel_spec(σᵤ²(zmat), ...)
 ```
 """
 function σᵤ²(arg::Vararg)
     return (:σᵤ², collect(:($(arg))))
 end

 """
 σᵤ²(arg::Vector)

An argument in `sfmodel_init()`. Specify the initial values for coefficients in
the `σᵤ²` function.

See the help on `sfmodel_init()` for more information.

# Examples
```julia-repl
sfmodel_init( σᵤ²(0.1, 0.2, 0.5), ...)

b0 = ones(3)*0.1
sfmodel_init( σᵤ²(b0), ...)
```
"""
 function σᵤ²(arg::Vector)
     return (:σᵤ², arg)
 end

  # -------------------

  """
  @sigma_u_2(arg::Vararg)

  alias of @σᵤ². See help on @σᵤ².
  """
 macro sigma_u_2(arg::Vararg)
     return (:σᵤ², collect(:($(arg))))
 end

 """
 sigma_u_2(arg::Vararg)

 alias of σᵤ². See help on σᵤ².
 """
 function sigma_u_2(arg::Vararg)
     return (:σᵤ², collect(:($(arg))))
 end

 """
 sigma_u_2(arg::Vector)

 alias of σᵤ². See help on σᵤ².
 """
 function sigma_u_2(arg::Vector)
     return (:σᵤ², arg)
 end

  # -------------------

  """
  @σᵥ²(arg::Vararg)

  An argument in `sfmodel_sepc()`. Specify the variables in the `σᵥ²`
  function using column names from a DataFrame.

  See the help on `sfmodel_spec()` for more information.

  # Examples
  ```julia-repl
  julia> df
  100×5 DataFrame
  │ Row │ yvar  │ xvar1     │ xvar2     │ zvar      │ _cons   │
  │     │ Int64 │ Float64   │ Float64   │ Float64   │ Float64 │
  ├─────┼───────┼───────────┼───────────┼───────────┼─────────┤
  │ 1   │ 1     │ 0.0306449 │ 0.452148  │ 0.808817  │ 1.0     │
  │ 2   │ 2     │ 0.460691  │ 0.296092  │ 0.454545  │ 1.0     │
  │ 3   │ 3     │ 0.897503  │ 0.376972  │ 0.907454  │ 1.0     │
  │ 4   │ 4     │ 0.682894  │ 0.776861  │ 0.161721  │ 1.0     │
  ⋮
  │ 96  │ 96    │ 0.329647  │ 0.0914057 │ 0.825032  │ 1.0     │
  │ 97  │ 97    │ 0.0781165 │ 0.338999  │ 0.761652  │ 1.0     │
  │ 98  │ 98    │ 0.41394   │ 0.0063118 │ 0.295372  │ 1.0     │
  │ 99  │ 99    │ 0.516381  │ 0.285415  │ 1.91995   │ 1.0     │
  │ 100 │ 100   │ 0.944     │ 0.702226  │ -0.539848 │ 1.0     │

  sfmodel_spec( @σᵥ²(zvar, _cons), ...)
  ```
  """
 macro σᵥ²(arg::Vararg)
     return (:σᵥ², collect(:($(arg))))
 end

 """
 σᵥ²(arg::Vararg)

 An argument in `sfmodel_sepc()`. Specify the name of the maxrix
 used as the data of the `σᵥ²` function.

 See the help on `sfmodel_spec()` for more information.

 # Examples
 ```julia-repl
 julia> sigv2mat
 100×2 Matrix{Float64}:
  -0.943133   1.0
  -0.897392   1.0
   0.585447   1.0
  -0.46106    1.0
  -0.54563    1.0
  -0.619428   1.0
   0.0575559  1.0
   ⋮
   0.0844192  1.0
  -1.3339     1.0
   1.29332    1.0
   0.691466   1.0
   0.422962   1.0
   0.374425   1.0

 sfmodel_spec(σᵥ²(sigv2mat), ...)
 ```
 """
 function σᵥ²(arg::Vararg)
     return (:σᵥ², collect(:($(arg))))
 end

 """
 σᵥ²(arg::Vector)

An argument in `sfmodel_init()`. Specify the initial values for coefficients in
the `σᵥ²` function.

See the help on `sfmodel_init()` for more information.

# Examples
```julia-repl
sfmodel_init( σᵥ²(0.1, 0.2, 0.5), ...)

b0 = ones(3)*0.1
sfmodel_init( σᵥ²(b0), ...)
```
"""
 function σᵥ²(arg::Vector)
     return (:σᵥ², arg)
 end

  # -------------------

"""
  @sigma_v_2(arg::Vararg)

 alias of @σᵥ². See help on @σᵥ².
"""
 macro sigma_v_2(arg::Vararg)
     return (:σᵥ², collect(:($(arg))))
 end

 """
  sigma_v_2(arg::Vararg)

 alias of σᵥ². See help on σᵥ².
 """
 function sigma_v_2(arg::Vararg)
     return (:σᵥ², collect(:($(arg))))
 end

 """
 sigma_v_2(arg::Vector)

alias of σᵥ². See help on σᵥ².
"""
 function sigma_v_2(arg::Vector)
     return (:σᵥ², arg)
 end

 # -------------------

 """
 @σₐ²(arg::Vararg)

 An argument in `sfmodel_sepc()`. Specify the variables in the `σₐ²`
 function using column names from a DataFrame.

 See the help on `sfmodel_spec()` for more information.

 # Examples
 ```julia-repl
 julia> df
 100×5 DataFrame
 │ Row │ yvar  │ xvar1     │ xvar2     │ zvar      │ _cons   │
 │     │ Int64 │ Float64   │ Float64   │ Float64   │ Float64 │
 ├─────┼───────┼───────────┼───────────┼───────────┼─────────┤
 │ 1   │ 1     │ 0.0306449 │ 0.452148  │ 0.808817  │ 1.0     │
 │ 2   │ 2     │ 0.460691  │ 0.296092  │ 0.454545  │ 1.0     │
 │ 3   │ 3     │ 0.897503  │ 0.376972  │ 0.907454  │ 1.0     │
 │ 4   │ 4     │ 0.682894  │ 0.776861  │ 0.161721  │ 1.0     │
 ⋮
 │ 96  │ 96    │ 0.329647  │ 0.0914057 │ 0.825032  │ 1.0     │
 │ 97  │ 97    │ 0.0781165 │ 0.338999  │ 0.761652  │ 1.0     │
 │ 98  │ 98    │ 0.41394   │ 0.0063118 │ 0.295372  │ 1.0     │
 │ 99  │ 99    │ 0.516381  │ 0.285415  │ 1.91995   │ 1.0     │
 │ 100 │ 100   │ 0.944     │ 0.702226  │ -0.539848 │ 1.0     │

 sfmodel_spec( @σₐ²(_cons), ...)
 ```
 """
macro σₐ²(arg::Vararg)
    return (:σₐ², collect(:($(arg))))
end

"""
σₐ²(arg::Vararg)

An argument in `sfmodel_sepc()`. Specify the name of the maxrix
used as the data of the `σₐ²` function.

See the help on `sfmodel_spec()` for more information.

# Examples
```julia-repl
julia> amat
100×2 Matrix{Float64}:
 -0.943133   1.0
 -0.897392   1.0
  0.585447   1.0
 -0.46106    1.0
 -0.54563    1.0
 -0.619428   1.0
  0.0575559  1.0
  ⋮
  0.0844192  1.0
 -1.3339     1.0
  1.29332    1.0
  0.691466   1.0
  0.422962   1.0
  0.374425   1.0

sfmodel_spec(σₐ²(amat), ...)
```
"""
function σₐ²(arg::Vararg)
    return (:σₐ², collect(:($(arg))))
end

"""
σₐ²(arg::Vector)

An argument in `sfmodel_init()`. Specify the initial values for coefficients in
the `σₐ²` function.

See the help on `sfmodel_init()` for more information.

# Examples
```julia-repl
sfmodel_init( σₐ²(0.1, 0.2, 0.5), ...)

b0 = ones(3)*0.1
sfmodel_init( σₐ²(b0), ...)
```
"""
function σₐ²(arg::Vector)
    return (:σₐ², arg)
end

 # -------------------

 """
 @sigma_a_2(arg::Vararg)

 alias of @σₐ². See help on @σₐ².
 """
macro sigma_a_2(arg::Vararg)
    return (:σₐ², collect(:($(arg))))
end

"""
sigma_a_2(arg::Vararg)

alias of σₐ². See help on σₐ².
"""
function sigma_a_2(arg::Vararg)
    return (:σₐ², collect(:($(arg))))
end

"""
sigma_a_2(arg::Vector)

alias of σₐ². See help on σₐ².
"""
function sigma_a_2(arg::Vector)
    return (:σₐ², arg)
end

  # -------------------

#=

 macro λ(arg::Vararg)
     return (:λ, collect(:($(arg))))
 end

 function λ(arg::Vararg)
     return (:λ, collect(:($(arg))))
 end

 function λ(arg::Vector)
     return (:λ, arg)
 end

 =#

  # -------------------

#=

 macro lambda(arg::Vararg)
     return (:λ, collect(:($(arg))))
 end

 function lambda(arg::Vararg)
     return (:λ, collect(:($(arg))))
 end

 function lambda(arg::Vector)
     return (:λ, arg)
 end

=#

# -------------------

#=

macro τ(arg::Vararg)
    return (:τ, collect(:($(arg))))
end

function τ(arg::Vararg)
    return (:τ, collect(:($(arg))))
end

function τ(arg::Vector)
    return (:τ, arg)
end

=#

 # -------------------

#=

macro tau(arg::Vararg)
    return (:τ, collect(:($(arg))))
end

function tau(arg::Vararg)
    return (:τ, collect(:($(arg))))
end

function tau(arg::Vector)
    return (:τ, arg)
end

=#

 # -------------------

 """
 @hscale(arg::Vararg)

 An argument in `sfmodel_sepc()`. Specify the variables in the `hscale`
 function using column names from a DataFrame.

 See the help on `sfmodel_spec()` for more information.

 # Examples
 ```julia-repl
 julia> df
 100×5 DataFrame
 │ Row │ yvar  │ xvar1     │ xvar2     │ zvar      │ _cons   │
 │     │ Int64 │ Float64   │ Float64   │ Float64   │ Float64 │
 ├─────┼───────┼───────────┼───────────┼───────────┼─────────┤
 │ 1   │ 1     │ 0.0306449 │ 0.452148  │ 0.808817  │ 1.0     │
 │ 2   │ 2     │ 0.460691  │ 0.296092  │ 0.454545  │ 1.0     │
 │ 3   │ 3     │ 0.897503  │ 0.376972  │ 0.907454  │ 1.0     │
 │ 4   │ 4     │ 0.682894  │ 0.776861  │ 0.161721  │ 1.0     │
 ⋮
 │ 96  │ 96    │ 0.329647  │ 0.0914057 │ 0.825032  │ 1.0     │
 │ 97  │ 97    │ 0.0781165 │ 0.338999  │ 0.761652  │ 1.0     │
 │ 98  │ 98    │ 0.41394   │ 0.0063118 │ 0.295372  │ 1.0     │
 │ 99  │ 99    │ 0.516381  │ 0.285415  │ 1.91995   │ 1.0     │
 │ 100 │ 100   │ 0.944     │ 0.702226  │ -0.539848 │ 1.0     │

 sfmodel_spec( @hscale(zvar, _cons), ...)
 ```
 """
macro hscale(arg::Vararg)
    return (:hscale, collect(:($(arg))))
end

"""
hscale(arg::Vararg)

An argument in `sfmodel_sepc()`. Specify the name of the maxrix
used as the data of the `hscale` function.

See the help on `sfmodel_spec()` for more information.

# Examples
```julia-repl
julia> hmat
100×2 Matrix{Float64}:
 -0.943133   1.0
 -0.897392   1.0
  0.585447   1.0
 -0.46106    1.0
 -0.54563    1.0
 -0.619428   1.0
  0.0575559  1.0
  ⋮
  0.0844192  1.0
 -1.3339     1.0
  1.29332    1.0
  0.691466   1.0
  0.422962   1.0
  0.374425   1.0

sfmodel_spec(hscale(hmat), ...)
```
"""
function hscale(arg::Vararg)
    return (:hscale, collect(:($(arg))))
end

"""
hscale(arg::Vector)

An argument in `sfmodel_init()`. Specify the initial values for coefficients in
the `hscale` function.

See the help on `sfmodel_init()` for more information.

# Examples
```julia-repl
sfmodel_init( hscale(0.1, 0.2, 0.5), ...)

b0 = ones(3)*0.1
sfmodel_init( hscale(b0), ...)
```
"""
function hscale(arg::Vector)
    return (:hscale, arg)
end

 # ---- time decay model ---------------

 """
 @gamma(arg::Vararg)

 An argument in `sfmodel_sepc()`. Specify the variables in the `gamma`
 function using column names from a DataFrame.

 See the help on `sfmodel_spec()` for more information.

 # Examples
 ```julia-repl
 julia> df
 100×5 DataFrame
 │ Row │ yvar  │ xvar1     │ xvar2     │ zvar      │ _cons   │
 │     │ Int64 │ Float64   │ Float64   │ Float64   │ Float64 │
 ├─────┼───────┼───────────┼───────────┼───────────┼─────────┤
 │ 1   │ 1     │ 0.0306449 │ 0.452148  │ 0.808817  │ 1.0     │
 │ 2   │ 2     │ 0.460691  │ 0.296092  │ 0.454545  │ 1.0     │
 │ 3   │ 3     │ 0.897503  │ 0.376972  │ 0.907454  │ 1.0     │
 │ 4   │ 4     │ 0.682894  │ 0.776861  │ 0.161721  │ 1.0     │
 ⋮
 │ 96  │ 96    │ 0.329647  │ 0.0914057 │ 0.825032  │ 1.0     │
 │ 97  │ 97    │ 0.0781165 │ 0.338999  │ 0.761652  │ 1.0     │
 │ 98  │ 98    │ 0.41394   │ 0.0063118 │ 0.295372  │ 1.0     │
 │ 99  │ 99    │ 0.516381  │ 0.285415  │ 1.91995   │ 1.0     │
 │ 100 │ 100   │ 0.944     │ 0.702226  │ -0.539848 │ 1.0     │

 sfmodel_spec( @gamma(zvar, _cons), ...)
 ```
 """
macro gamma(arg::Vararg)
    return (:gamma, collect(:($(arg))))
end

"""
gamma(arg::Vararg)

An argument in `sfmodel_sepc()`. Specify the name of the maxrix
used as the data of the `gamma` function.

See the help on `sfmodel_spec()` for more information.

# Examples
```julia-repl
julia> gmat
100×2 Matrix{Float64}:
 -0.943133   1.0
 -0.897392   1.0
  0.585447   1.0
 -0.46106    1.0
 -0.54563    1.0
 -0.619428   1.0
  0.0575559  1.0
  ⋮
  0.0844192  1.0
 -1.3339     1.0
  1.29332    1.0
  0.691466   1.0
  0.422962   1.0
  0.374425   1.0

sfmodel_spec(gamma(gmat), ...)
```
"""
function gamma(arg::Vararg)
    return (:gamma, collect(:($(arg))))
end

"""
gamma(arg::Vector)

An argument in `sfmodel_init()`. Specify the initial values for coefficients in
the `gamma` function.

See the help on `sfmodel_init()` for more information.

# Examples
```julia-repl
sfmodel_init( gamma(0.1, 0.2, 0.5), ...)

b0 = ones(3)*0.1
sfmodel_init( gamma(b0), ...)
```
"""
function gamma(arg::Vector)
    return (:gamma, arg)
end

 # ------ panel model -------------

  macro sfpanel(arg::Vararg)

    (length(arg)==1)||throw("@sfpanel in sfmodel_spec should have one single string.")

    if !(arg[1] ∈ (:TFE_WH2010, :TFE_CSW2014, :TFE_G2005, :TRE, :TimeDecay, :Kumbhakar1990 ))
         throw("The keyword of @sfpanel in sfmodel_spec is specified incorrectly.")
    end

    if arg[1] == :TFE_G2005
        throw("The TEF_G2005 panel model is essentially a cross-sectional model with dummies of individuals. Please use non-panel specifications and add individual dummies to @frontier().")
    end

    return (:panel, collect(:($(arg))))
 end

 function sfpanel(arg::Vararg)

    (length(arg)==1) || throw("`sfpanel()` in sfmodel_spec should have one single string.")

    if !(arg[1] ∈ (TFE_WH2010, TFE_CSW2014, TFE_G2005, TRE, TimeDecay, Kumbhakar1990 ))
         throw("The keyword of `sfpanel` in sfmodel_spec is specified incorrectly.")
    end

    if arg[1] == TFE_G2005
        throw("The TEF_G2005 panel model is essentially a cross-sectional model with dummies of individuals. Please use non-panel specifications and add individual dummies to @frontier().")
    end

    return (:panel, [Symbol(arg[1])])
 end

 """
 @timevar(arg::Vararg)

 An argument in `sfmodel_sepc()` for panel data models. Specify the column
 name of a DataFrame that contain the time information of the panel data.

 See the help on `sfmodel_spec()` for more information.

 # Examples
 ```julia-repl
 julia> df
 100×4 DataFrame
 │ Row │ year  │ firm  │ yvar  │ xvar1     │
 │     │ Int64 │ Int64 │ Int64 │ Float64   │
 ├─────┼───────┼───────┼───────┼───────────┤
 │ 1   │ 2019  │ 1     │ 1     │ 0.77645   │
 │ 2   │ 2020  │ 1     │ 2     │ 0.0782388 │
 │ 3   │ 2021  │ 1     │ 3     │ 0.222884  │
 │ 4   │ 2022  │ 1     │ 4     │ 0.762864  │
 ⋮
 │ 96  │ 2022  │ 24    │ 96    │ 0.590184  │
 │ 97  │ 2019  │ 25    │ 97    │ 0.364425  │
 │ 98  │ 2020  │ 25    │ 98    │ 0.639463  │
 │ 99  │ 2021  │ 25    │ 99    │ 0.500526  │
 │ 100 │ 2022  │ 25    │ 100   │ 0.239137  │

 sfmodel_spec( @timvar(year), ...)
 ```
 """
macro timevar(arg::Vararg)
    return (:timevar, collect(:($(arg))))
end

"""
timevar(arg::Vararg)

An argument in `sfmodel_sepc()` for panel data models. Specify the name of the
matrix containing the time information of the panel data.

See the help on `sfmodel_spec()` for more information.

# Examples
```julia-repl
julia> timemat
100-element Vector{Int64}:
 2019
 2020
 2021
 2022
 2019
 2020
 2021
    ⋮
 2021
 2022
 2019
 2020
 2021
 2022

sfmodel_spec( timvar(timemat), ...)
```
"""
function timevar(arg::Vararg)
    return (:timevar, collect(:($(arg))))
end

"""
@idvar(arg::Vararg)

An argument in `sfmodel_sepc()` for panel data models. Specify the column
name of a DataFrame that contain the individual's id information of the panel data.

See the help on `sfmodel_spec()` for more information.

# Examples
```julia-repl
julia> df
100×4 DataFrame
│ Row │ year  │ firm  │ yvar  │ xvar1     │
│     │ Int64 │ Int64 │ Int64 │ Float64   │
├─────┼───────┼───────┼───────┼───────────┤
│ 1   │ 2019  │ 1     │ 1     │ 0.77645   │
│ 2   │ 2020  │ 1     │ 2     │ 0.0782388 │
│ 3   │ 2021  │ 1     │ 3     │ 0.222884  │
│ 4   │ 2022  │ 1     │ 4     │ 0.762864  │
⋮
│ 96  │ 2022  │ 24    │ 96    │ 0.590184  │
│ 97  │ 2019  │ 25    │ 97    │ 0.364425  │
│ 98  │ 2020  │ 25    │ 98    │ 0.639463  │
│ 99  │ 2021  │ 25    │ 99    │ 0.500526  │
│ 100 │ 2022  │ 25    │ 100   │ 0.239137  │

sfmodel_spec( @idvar(firm), ...)
```
"""
macro idvar(arg::Vararg)
    return (:idvar, collect(:($(arg))))
end

"""
idvar(arg::Vararg)

An argument in `sfmodel_sepc()` for panel data models. Specify the name of the
matrix containing the individual's id information of the panel data.

See the help on `sfmodel_spec()` for more information.

# Examples
```julia-repl
julia> idmat
100-element Vector{Int64}:
  1
  1
  1
  1
  2
  2
  2
  ⋮
 24
 24
 25
 25
 25
 25

sfmodel_spec( idvar(idmat), ...)
```
"""
function idvar(arg::Vararg)
    return (:idvar, collect(:($(arg))))
end

function misc(arg::Any)
    return (:misc, arg)
end

     #* ------  variables by vector type, not in use -----
 #=
 macro xvar(arg::Vararg)
     return (:xvar, collect(:($(arg))))
 end

  # -------------------

 macro zvar(arg::Vararg)
     return (:zvar, collect(:($(arg))))
 end

  # -------------------

 macro etavar(arg::Vararg)
     return (:etavar, collect(:($(arg))))
 end
 =#

 function all_init(arg::Vararg)
     return (:all_init, collect(:($(arg))))
 end

"""
all_init(arg::Vector)

An argument in `sfmodel_init()`. Specify the initial values for all the
parameters in the model.

See the help on `sfmodel_init()` for more information.

# Examples
```julia-repl
sfmodel_init( all_init(0.1, 0.2, 0.5, -0.1, -0.1, -0.2), ...)

b0 = ones(6)*0.1
sfmodel_init( all_init(b0), ...)
```
"""
 function all_init(arg::Vector)
     return (:all_init, arg)
 end

 #? ----  functions for sfmodel_opt --------

 #= macro warmstart_solver(arg=nothing)
     return (:warmstart_solver, arg)
 end =#

"""
warmstart_solver(arg)

An argument in `sfmodel_opt()`. Specify the algorithm used in the first-stage ("warmstart")
  optimization process.

The default is `NelderMead()`. Others include `SimulatedAnnealing()`, `SAMIN()`, `ParticleSwarm()`,
  `ConjugateGradient()`, `GradientDescent()`, `BFGS()`, `LBFGS()`,
  `Newton()`, `NewtonTrustRegion()`, and `IPNewton()`. See
  http://julianlsolvers.github.io/Optim.jl/stable/ for details.
  Non-gradient based algorithms are recommended for the warmstart solver.

See the help on `sfmodel_opt()` for more information.

# Examples
```julia-repl
sfmodel_opt( warmstart_solver(NelderMead()), ...)
```
"""
 function warmstart_solver(arg=nothing)
     return (:warmstart_solver, arg)
 end

  # -------------------

 #= macro warmstart_maxIT(arg=nothing)
     return (:warmstart_maxIT, arg)
 end =#

 """
warmstart_maxIT(arg)

An argument in `sfmodel_opt()`. Specify the iteration limit for the warmstart. Default
is 100.

See the help on `sfmodel_opt()` for more information.

# Examples
```julia-repl
sfmodel_opt( warmstart_maxIT(400), ...)
```
"""
 function warmstart_maxIT(arg=nothing)
     return (:warmstart_maxIT, arg)
 end

  # -------------------

 #= macro main_solver(arg=nothing)
     return (:main_solver, arg)
 end =#

 """
main_solver(arg)

An argument in `sfmodel_opt()`. Specify the algorithm used in the 2nd-stage ("main")
  optimization process.

  The default is `Newton()`. Others include `SimulatedAnnealing()`, `SAMIN()`, `ParticleSwarm()`,
  `ConjugateGradient()`, `GradientDescent()`, `BFGS()`, `LBFGS()`,
  `NewtonTrustRegion()`, and `IPNewton()`. See
  http://julianlsolvers.github.io/Optim.jl/stable/ for details.

See the help on `sfmodel_opt()` for more information.

# Examples
```julia-repl
sfmodel_opt( main_solver(Newton()), ...)
```
"""
 function main_solver(arg=nothing)
     return (:main_solver, arg)
 end

  # -------------------

 #= macro main_maxIT(arg=nothing)
     return (:main_maxIT, arg)
 end =#

 """
main_maxIT(arg)

 An argument in `sfmodel_opt()`. Specify the iteration limit for the main solver. Default
 is 2000.

 See the help on `sfmodel_opt()` for more information.

 # Examples
 ```julia-repl
 sfmodel_opt( main_maxIT(2500), ...)
 ```
 """
 function main_maxIT(arg::Any=nothing)
     return (:main_maxIT, arg)
 end

  # -------------------

 #= macro tolerance(arg=nothing)
     return (:tolerance, arg)
 end =#

 """
tolerance(arg::Float64)

An argument in `sfmodel_opt()`. Specify the convergence criterion ("tolerance") based on the
absolute value of gradients. Default is 1.0e-8. For non-gradient algorithms,
it controls the main convergence tolerance, which is solver specific.
See `Optim`'s `g_tol` option for more information.

Also see the help on `sfmodel_opt()` for more information.

# Examples
```julia-repl
sfmodel_opt( tolerance(1.0e-6), ...)
```
"""
 function tolerance(arg::Float64=nothing)
     return (:tolerance, arg)
 end

   # --------------------------
 """
verbose(arg::Bool)

An argument in `sfmodel_opt()`. Specify whether to print on screen (`true`,
the default) the information of the model and the optimization results.

See the help on `sfmodel_opt()` for more information.

# Examples
```julia-repl
sfmodel_opt( verbose(false), ...)
```
"""
 function verbose(arg::Bool=false)
     return (:verbose, arg)
 end

 """
 banner(arg::Bool)

 An argument in `sfmodel_opt()`. Specify whether to print on screen (`true`,
 the default) a banner to serve as a visual indicator of the start of the
 estimation.

 See the help on `sfmodel_opt()` for more information.

 # Examples
 ```julia-repl
 sfmodel_opt( banner(false), ...)
 ```
 """
 function banner(arg::Bool=true)
    return (:banner, arg)
 end

 """
 ineff_index(arg::Bool)

 An argument in `sfmodel_opt()`. Specify whether (`true`, the default) to compute the Jondrow et al. (1982)
 inefficiency index and the Battese and Coelli (1988) efficiency index.

 See the help on `sfmodel_opt()` for more information.

 # Examples
 ```julia-repl
 sfmodel_opt( ineff_index(false), ...)
 ```
 """
 function ineff_index(arg::Bool=true)
    return (:ineff_index, arg)
 end

"""
marginal(arg::Bool)

 An argument in `sfmodel_opt()`. Specify whether (`true`, the default) to
 compute the marginal effects of the exogenous determinants of inefficiency (if any).

 See the help on `sfmodel_opt()` for more information.

 # Examples
 ```julia-repl
 sfmodel_opt( marginal(false), ...)
 ```
"""
 function marginal(arg::Bool=true)
    return (:marginal, arg)
 end

  # -------------------------

  """
  table_format(arg)

   An argument in `sfmodel_opt()`. Specify the format to print the coefficient
   tables on the screen: `text` (default), `html`, or `latex`. A wrapper of `PrettyTables.jl`'s
   `backend` option.

   See the help on `sfmodel_opt()` for more information.

   # Examples
   ```julia-repl
   sfmodel_opt( table_format(html), ...)
   ```
  """
 function table_format(arg=nothing)

    if !(arg ∈ (text, html, latex ))
        throw("The keyword of `table_format` in `sfmodel_opt()` is specified incorrectly. Only allow `text`, `html`, or `latex`. Got `$(arg)` instead.")
   end

    return (:table_format, Symbol(arg))
end

  #? ------ functions for sfmodel_fit ----------

 # macro dataframe(arg)
 #     return arg
 # end

 """
 useData(D::DataFrame)

  An argument in `sfmodel_fit()`. Specify the name of the DataFrame that
  contains the estimation data.

  See the help on `sfmodel_fit()` for more information.

# Examples
```julia-repl
  julia> mydf
  100×4 DataFrame
  │ Row │ year  │ firm  │ yvar  │ xvar1     │
  │     │ Int64 │ Int64 │ Int64 │ Float64   │
  ├─────┼───────┼───────┼───────┼───────────┤
  │ 1   │ 2019  │ 1     │ 1     │ 0.77645   │
  │ 2   │ 2020  │ 1     │ 2     │ 0.0782388 │
  │ 3   │ 2021  │ 1     │ 3     │ 0.222884  │
  │ 4   │ 2022  │ 1     │ 4     │ 0.762864  │
  ⋮
  │ 96  │ 2022  │ 24    │ 96    │ 0.590184  │
  │ 97  │ 2019  │ 25    │ 97    │ 0.364425  │
  │ 98  │ 2020  │ 25    │ 98    │ 0.639463  │
  │ 99  │ 2021  │ 25    │ 99    │ 0.500526  │
  │ 100 │ 2022  │ 25    │ 100   │ 0.239137  │

  sfmodel_fit(useData(mydf), ...)
```
"""
 function useData(D::DataFrame)  # doesn't work using macro (perhaps because of data), so...
     return D
 end

 #? ---- macros for sf_predict -------------

 """
@eq(arg)

An argument in `sfmodel_predict()`. Specify the name of the function to be
predicted.

See the help on `sfmodel_predict()` for more information.

# Examples
```julia-repl
sfmodel_predict( @eq(frontier), ...)
```
"""
 macro eq(arg)
    return [:($(arg))]
end


text:
@dict_readable struct WorkspaceFoldersChangeEvent
    added::Vector{WorkspaceFolder}
    removed::Vector{WorkspaceFolder}
end

@dict_readable struct DidChangeWorkspaceFoldersParams
    event::WorkspaceFoldersChangeEvent
end

@dict_readable struct DidChangeConfigurationParams
    settings::Any
end

struct ConfigurationItem <: Outbound
    scopeUri::Union{DocumentUri,Missing}
    section::Union{String,Missing}
end

struct ConfigurationParams <: Outbound
    items::Vector{ConfigurationItem}
end

##############################################################################
# File watching
const FileChangeType = Int
const FileChangeTypes = (Created = 1,
                         Changed = 2,
                         Deleted = 3)

const WatchKind = Int
const WatchKinds = (Create = 1,
                    Change = 2,
                    Delete = 4)

@dict_readable struct FileEvent
    uri::String
    type::FileChangeType
end

@dict_readable struct DidChangeWatchedFilesParams
    changes::Vector{FileEvent}
end

struct FileSystemWatcher <: Outbound
    globPattern::String
    kind::Union{WatchKind,Missing}
end

struct DidChangeWatchedFilesRegistrationOptions <: Outbound
    watchers::Vector{FileSystemWatcher}
end

##############################################################################
# Registration

struct Registration <: Outbound
    id::String
    method::String
    registerOptions::Union{Any,Missing}
end

struct RegistrationParams <: Outbound
    registrations::Vector{Registration}
end

struct TextDocumentRegistrationOptions <: Outbound
    documentSelector::Union{DocumentSelector,Nothing}
end

struct TextDocumentChangeRegistrationOptions <: Outbound
    documentSelector::DocumentSelector
    syncKind::Int
end

struct TextDocumentSaveRegistrationOptions <: Outbound
    documentSelector::DocumentSelector
    includeText::Union{Bool,Missing}
end

struct Unregistration <: Outbound
    id::String
    method::String
end

struct UnregistrationParams <: Outbound
    unregistrations::Vector{Unregistration}
end

struct StaticRegistrationOptions <: Outbound
    id::Union{String,Missing}
end

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


text:
include("TestBaselines/JointPlot.jl")
include("TestBaselines/Raw.jl")
include("TestBaselines/Standardised.jl")
include("TestBaselines/Filter.jl")
include("TestBaselines/Transform.jl")
include("TestBaselines/IntervalScale.jl")
include("TestBaselines/Heuristic.jl")
include("TestBaselines/bPCA.jl")
include("TestBaselines/RobustIntervalScale.jl")


text:
################################################################################

"""

  mindGraphics(modelHMM::Dict{S, Tuple{Array{T, 1}, Array{Array{U, 1}, 1}}}, shArgs::Dict;
  labels = nothing) where S <: String where T <: Int64 where U <: Float64

# Description
Control graphic module of *MindReader*. Verifies arguments, builds heatmap, write to file and renders image.

# Arguments
`modelHMM` model to render by graphics module.

`shArgs` shell arguments.

if available `labels` are available, concatenate to rendering.

"""
function mindGraphics(modelHMM::Dict{S, Tuple{Array{T, 1}, Array{Array{U, 1}, 1}}}, shArgs::Dict; labels = nothing) where S <: String where T <: Int64 where U <: Float64

  # check arguments | assign
  κ = [:svg, :csv]
  for (ζ, ξ) ∈ zip([Symbol("out", ι) for ι ∈ κ], κ)
    τ = shChecker(shArgs, string(ζ), string(ξ))
    @eval $ζ = $τ
  end

  @info "Plotting..."
  # create array to plot
  toHeat, keyString = collectState(modelHMM)

  if !isnothing(labels)
    # concatenate annotations
    toHeat = [toHeat; labels' .+ 1]
  end

  # write
  writedlm( string(shArgs["outdir"], outcsv), toHeat, "," )

  @info "Rendering..."
  renderGraphics( string(shArgs["outdir"], outsvg), toHeat )

end

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

"wrapper for model heatmap plotting"
function renderGraphics(filename::S, toHeat::Array{T, 2}) where S <: String where T <: Number

  #  TODO: add channel labels by passing vector to heatmap function
  # plot layout
  plotFig = Figure()
  heatplot = plotFig[1, 1] = Axis(plotFig, title = "States Heatmap")

  # axis labels
  heatplot.xlabel = "Recording Time"
  heatplot.ylabel = "Channels"
  # heatplot.yticks = 1:22

  # heatmap plot & color range
  hm = heatmap!(heatplot, toHeat')
  hm.colorrange = (0, 5)

  # color bar
  cbar = plotFig[2, 1] = Colorbar(plotFig, hm, label = "HMM states")
  cbar.vertical = false
  cbar.height = 10
  cbar.width = Relative(2 / 3)
  cbar.ticks = 1:1:5

  # save rendering
  save( filename, plotFig )

end

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

"iterate over known electrodes and collect states from Hidden Markov model"
function collectState(modelHMM::Dict{S, Tuple{Array{T, 1}, Array{Array{U, 1}, 1}}}) where S <: String where T <: Int64 where U <: Float64

  keyString = modelHMM[convert.(String, keys(modelHMM))[1]][1]
  toHeat = zeros(length(modelHMM), length(keyString))
  ψ = size(toHeat, 1)
  for ε ∈ elecID
    if haskey(modelHMM, ε)
      toHeat[ψ, :] = modelHMM[ε][1]
      ψ -= 1
    else
      @debug ε
    end
  end
  return toHeat, keyString
end

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

"check whether shell argument was explicitly declared and assigned"
function shChecker(shArgs::Dict, κ::S, ζ::S) where S <: String
  if shArgs[κ] != nothing
    return shArgs[κ]
  else
    return replace(shArgs["file"], "edf" => ζ)
  end
end

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

max_stars_repo_name: DanielRivasMD/MindReader

max_issues_repo_name: DanielRivasMD/MindReader

max_forks_repo_name: DanielRivasMD/MindReader


text:
using jInv.Mesh
using jInv.LinearSolvers
using Multigrid.ParallelJuliaSolver
using Multigrid.DomainDecomposition
using Multigrid
using KrylovMethods

using SparseArrays
using LinearAlgebra
using Test
using jInv.Mesh;

include("DDPoissonFuncs.jl");

println("************************************************* Dirichlet DomainDecomposition for Poisson 2D ******************************************************");

domain = [0.0, 1.0, 0.0, 1.0];
n = [32,32];
Mr = getRegularMesh(domain,n)
G = getNodalGradientMatrix(Mr);
Ar = G'*G;
Ar = Ar + 1e-5*opnorm(Ar,1)*sparse(1.0I,size(Ar,2),size(Ar,2));

NumCells = [8,8];
overlap = [1,1];
n_tup = tuple(n...);
dim = length(n);

numCores 	= 2;
maxIter     = 10;
relativeTol = 1e-6;

q = randn(prod(n.+1));
x = zeros(size(q));
Ainv = getParallelJuliaSolver(Float64,Int64,numCores=numCores,backend=1);
# Ainv = getJuliaSolver();
DDparam = getDomainDecompositionParam(Mr,NumCells,overlap,getNodalIndicesOfCell,Ainv);

(~,DDparam) = solveLinearSystem!(Ar,q,x,DDparam)

println(norm(Ar'*x - q)/norm(q))

println("************************************************* Neumann DomainDecomposition for Poisson 2D ******************************************************");

Ctor = DomainDecompositionOperatorConstructor(Mr,getSubParams,getLap,getDirichletMassNodalMesh);

setupDDSerial(Ctor,DDparam);

Prec = r->solveDDSerial(Ar,r,zeros(eltype(q),size(q)),DDparam,1,0)[1];
x = zeros(size(q));

x, flag,rnorm,iter = KrylovMethods.fgmres(getAfun(Ar,zeros(eltype(q),size(q)),4),q,10,tol = 1e-10,maxIter = 2,M = Prec, x = x,out=2,flexible=true);

# x = zeros(size(q));

# x, flag,rnorm,iter = KrylovMethods.bicgstb(getAfun(Ar,zeros(eltype(q),size(q)),4), q, tol = 1e-10,maxIter = 25,M1 = Prec, x = x,out=2)

# println("****************************** DomainDecomposition for Poisson 3D ******************************")

# domain = [0.0, 1.0, 0.0, 1.0,0.0,1.0];
# n = [64,64,32];
# Mr = getRegularMesh(domain,n)
# G = getNodalGradientMatrix(Mr);
# Ar = G'*G;
# Ar = Ar + 1e-4*opnorm(Ar,1)*sparse(1.0I,size(Ar,2),size(Ar,2));
# N = size(Ar,2);

# q = Ar*rand(N);
# q = q/norm(q);
# x = zeros(N);

# numCores 	= 2;
# maxIter     = 10;
# relativeTol = 1e-6;

# NumCells = [4,4,2];
# overlap = [2,2,2];

# Ainv = getParallelJuliaSolver(Float64,Int64,numCores=numCores,backend=1);
# # Ainv = getJuliaSolver();
# DDparam = getDomainDecompositionParam(Mr,NumCells,overlap,getNodalIndicesOfCell,Ainv);

# (~,DDparam) = solveLinearSystem!(Ar,q,x,DDparam)
# # x = solveLinearSystem!(Ar,q,x,DDparam)[1];

# println(norm(Ar'*x - q)/norm(q))


text:
using Compat
using Compat.Libdl

const LIBMBHFMM2D = string(Base.Filesystem.dirname(Base.source_path()), "/../bin/libmbhfmm2d")

include("BoxFMMType.jl")
include("BoxFMMUtil.jl")

# Load OpenMP library
Libdl.dlopen("libgomp", Libdl.RTLD_GLOBAL)

mutable struct MBHFMM2DParams

    # storage of parameters which
    # describe the FMM tree and call

    lambda::Float64
    iprec::Int32
    ifalltarg::Bool
    src::Array{Float64,2}
    targ::Array{Float64,2}
    ifcharge::Bool
    ifdipole::Bool
    ifquad::Bool
    ifoct::Bool
    charge::Array{Float64}
    dipstr::Array{Float64}
    dipvec::Array{Float64,2}
    quadstr::Array{Float64}
    quadvec::Array{Float64,2}
    octstr::Array{Float64}
    octvec::Array{Float64,2}
    exrad::Array{Float64,1}
end

function MBHFMM2DParams(lambda::Float64,src::Array{Float64,2},
                        targ::Array{Float64,2},ifcharge::Bool,
                        ifdipole::Bool,ifquad::Bool,
                        ifoct::Bool,charge::Array{Float64,1},
                        dipstr::Array{Float64,1},
                        dipvec::Array{Float64,2},
                        quadstr::Array{Float64,1},
                        quadvec::Array{Float64,2},
                        octstr::Array{Float64,1},
                        octvec::Array{Float64,2};
                        iprec::Int=3,
                        ifalltarg::Bool=false,
                        exrad::Array{Float64,1}
                        =Array{Float64}(undef,0))

    iprec = convert(Int32,iprec)

    # constructor for MBHFMM2DParams, with a couple
    # default values

    return MBHFMM2DParams(lambda,iprec,ifalltarg,
                          src,targ,ifcharge,ifdipole,
                          ifquad,ifoct,
                          charge,dipstr,dipvec,quadstr,
                          quadvec,octstr,octvec,exrad)

end

mutable struct MBHFMM2DStorage

    # storage for FMM itself

    tree::BoxTree2D
    ifalltarg::Array{Int32,1}
    sorted_pts::SortedPts2D
    chargesort::Array{Float64}
    dipstrsort::Array{Float64}
    dipvecsort::Array{Float64,2}
    quadstrsort::Array{Float64}
    quadvecsort::Array{Float64,2}
    octstrsort::Array{Float64}
    octvecsort::Array{Float64,2}
    ifcharge::Array{Int32,1}
    ifdipole::Array{Int32,1}
    ifquad::Array{Int32,1}
    ifoct::Array{Int32,1}
    isave::Array{Int32,1}
    dsave::Array{Float64,1}
    csave::Array{Complex{Float64},1}
end

function MBHFMM2DStorage_init(fmmpars::MBHFMM2DParams,
                              tree::BoxTree2D,
                              sorted_pts::SortedPts2D)

    # convert bools to fortran flags

    ifalltarg = boxfmm2d_booltoint32(fmmpars.ifalltarg)
    ifcharge = boxfmm2d_booltoint32(fmmpars.ifcharge)
    ifdipole = boxfmm2d_booltoint32(fmmpars.ifdipole)
    ifquad = boxfmm2d_booltoint32(fmmpars.ifquad)
    ifoct = boxfmm2d_booltoint32(fmmpars.ifoct)

    # if charge is used, sort it

    if fmmpars.ifcharge
        chargesort = fmmpars.charge[sorted_pts.isrcsort]
    else
        chargesort = zeros(Float64,1)
    end

    if fmmpars.ifdipole
        dipstrsort = fmmpars.dipstr[sorted_pts.isrcsort]
        dipvecsort = fmmpars.dipvec[:,sorted_pts.isrcsort]
    else
        dipstrsort = zeros(Float64,1)
        dipvecsort = zeros(Float64,2,1)
    end

    if fmmpars.ifquad
        quadstrsort = fmmpars.quadstr[sorted_pts.isrcsort]
        quadvecsort = fmmpars.quadvec[:,sorted_pts.isrcsort]
    else
        quadstrsort = zeros(Float64,1)
        quadvecsort = zeros(Float64,3,1)
    end

    if fmmpars.ifoct
        octstrsort = fmmpars.octstr[sorted_pts.isrcsort]
        octvecsort = fmmpars.octvec[:,sorted_pts.isrcsort]
    else
        octstrsort = zeros(Float64,1)
        octvecsort = zeros(Float64,4,1)
    end

    # initialize isave, dsave, csave to short arrays

    isave = Array{Int32}(undef, 1)
    dsave = Array{Float64}(undef, 1)
    csave = Array{Complex{Float64}}(undef, 1)

    return MBHFMM2DStorage(tree,ifalltarg,sorted_pts,
                           chargesort,dipstrsort,
                           dipvecsort,quadstrsort,
                           quadvecsort,octstrsort,
                           octvecsort,ifcharge,ifdipole,
                           ifquad,ifoct,isave,dsave,csave)

end

function mbhfmm2d_form(fmmpars::MBHFMM2DParams;maxnodes::Int=30,
                       maxboxes::Int=-1)
    # form tree
    tree, sorted_pts, ier = mbhfmm2d_tree(fmmpars, maxnodes=maxnodes, maxboxes=maxboxes,exrad=fmmpars.exrad)
    # do work
    mbhfmm2d_form(fmmpars, tree, sorted_pts)
end

function mbhfmm2d_tree(fmmpars::MBHFMM2DParams;maxnodes::Int=30,
                       maxboxes::Int=-1,
                       exrad::Array{Float64,1}
                       =Array{Float64}(undef,0))
    if maxboxes < 0
        return [], 1
    end
    tree, sorted_pts, ier = BoxTree2DMaxBoxesST(fmmpars.src, fmmpars.targ, maxboxes,
                                                maxnodes=maxnodes,
                                                ifverbose=false,
                                                ifalltarg=fmmpars.ifalltarg,
                                                exrad=exrad)
    return tree, sorted_pts, ier
end

function mbhfmm2d_form(fmmpars::MBHFMM2DParams,
                       tree::BoxTree2D,
                       sorted_pts::SortedPts2D)

    # initialize storage

    fmmstor = MBHFMM2DStorage_init(fmmpars,tree,sorted_pts)

    # query fmm for length of isave, dsave, csave

    lambda = fmmpars.lambda
    ier1 = zeros(Int32,1)
    iprec = fmmpars.iprec
    nlev = tree.nlev
    levelbox = tree.levelbox
    iparentbox = tree.iparentbox
    ichildbox = tree.ichildbox
    icolbox = tree.icolbox
    irowbox = tree.irowbox
    nboxes = tree.nboxes
    nblevel = tree.nblevel
    iboxlev = tree.iboxlev
    istartlev = tree.istartlev
    ifalltarg = fmmstor.ifalltarg
    localonoff = tree.localonoff
    zll = tree.zll
    blength = tree.blength
    ns = sorted_pts.ns
    srcsort = sorted_pts.srcsort
    isrcladder = sorted_pts.isrcladder
    ifcharge = fmmstor.ifcharge
    chargesort = fmmstor.chargesort
    ifdipole = fmmstor.ifdipole
    dipstrsort = fmmstor.dipstrsort
    dipvecsort = fmmstor.dipvecsort
    ifquad = fmmstor.ifquad
    quadstrsort = fmmstor.quadstrsort
    quadvecsort = fmmstor.quadvecsort
    ifoct = fmmstor.ifoct
    octstrsort = fmmstor.octstrsort
    octvecsort = fmmstor.octvecsort
    isave = fmmstor.isave
    dsave = fmmstor.dsave
    csave = fmmstor.csave

    lisave = [convert(Int32,-1)]
    ldsave = [convert(Int32,-1)]
    lcsave = [convert(Int32,-1)]

    ccall( (:mbhfmm2d_form_,LIBMBHFMM2D), Cvoid,
           (Ref{Float64},Ref{Int32},Ref{Int32},Ref{Int32},
            Ref{Int32},Ref{Int32},Ref{Int32},Ref{Int32},
            Ref{Int32},Ref{Int32},Ref{Int32},Ref{Int32},
            Ref{Int32},Ref{Int32},Ref{Int32},Ref{Float64},
            Ref{Float64},Ref{Int32},Ref{Float64},Ref{Int32},
            Ref{Int32},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Float64},
            Ref{Int32},Ref{Int32},Ref{Float64},Ref{Int32},
            Ref{Complex{Float64}},Ref{Int32}),
           lambda,ier1,iprec,nlev,levelbox,iparentbox,
           ichildbox,icolbox,irowbox,nboxes,nblevel,
           iboxlev,istartlev,ifalltarg,localonoff,zll,
           blength,ns,srcsort,isrcladder,ifcharge,chargesort,
           ifdipole,dipstrsort,dipvecsort,ifquad,
           quadstrsort,quadvecsort,ifoct,octstrsort,
           octvecsort,isave,lisave,dsave,ldsave,csave,
           lcsave)

    # allocate

    fmmstor.isave = Array{Int32}(undef, lisave[1])
    fmmstor.dsave = Array{Float64}(undef, ldsave[1])
    fmmstor.csave = Array{Complex{Float64}}(undef, lcsave[1])

    isave = fmmstor.isave
    dsave = fmmstor.dsave
    csave = fmmstor.csave

    # form fmm

    ccall( (:mbhfmm2d_form_,LIBMBHFMM2D), Cvoid,
           (Ref{Float64},Ref{Int32},Ref{Int32},Ref{Int32},
            Ref{Int32},Ref{Int32},Ref{Int32},Ref{Int32},
            Ref{Int32},Ref{Int32},Ref{Int32},Ref{Int32},
            Ref{Int32},Ref{Int32},Ref{Int32},Ref{Float64},
            Ref{Float64},Ref{Int32},Ref{Float64},Ref{Int32},
            Ref{Int32},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Float64},
            Ref{Int32},Ref{Int32},Ref{Float64},Ref{Int32},
            Ref{Complex{Float64}},Ref{Int32}),
           lambda,ier1,iprec,nlev,levelbox,iparentbox,
           ichildbox,icolbox,irowbox,nboxes,nblevel,
           iboxlev,istartlev,ifalltarg,localonoff,zll,
           blength,ns,srcsort,isrcladder,ifcharge,chargesort,
           ifdipole,dipstrsort,dipvecsort,ifquad,
           quadstrsort,quadvecsort,ifoct,octstrsort,
           octvecsort,isave,lisave,dsave,ldsave,csave,
           lcsave)

    return fmmstor

end

function mbhfmm2d_targ!(fmmpars::MBHFMM2DParams,
                        fmmstor::MBHFMM2DStorage,
                        targ::Array{Float64,2},
                        ifpot::Bool,pottarg::Array{Float64,1},
                        ifgrad::Bool,gradtarg::Array{Float64,2},
                        ifhess::Bool,hesstarg::Array{Float64,2},
                        ifignore_exrad::Bool=false)

    # determine if an exclusion radius should be used
    # if ifignore_exrad == true, then the exclusion
    # radius is ignored

    sorted_pts = fmmstor.sorted_pts
    ifexrad = sorted_pts.ifexrad

    if (!ifexrad || ifignore_exrad)
        mbhfmm2d_targ_all!(fmmpars,fmmstor,targ,ifpot,pottarg,
                           ifgrad,gradtarg,ifhess,hesstarg)
    else
        exradsort = sorted_pts.exradsort
        mbhfmm2d_targ_exrad!(fmmpars,fmmstor,exradsort,
                             targ,ifpot,pottarg,
                             ifgrad,gradtarg,ifhess,hesstarg)
    end

end

function mbhfmm2d_targ_all!(fmmpars::MBHFMM2DParams,
                            fmmstor::MBHFMM2DStorage,
                            targ::Array{Float64,2},
                         ifpot::Bool,pottarg::Array{Float64,1},
                         ifgrad::Bool,gradtarg::Array{Float64,2},
                         ifhess::Bool,hesstarg::Array{Float64,2})

    # unpack

    tree = fmmstor.tree
    sorted_pts = fmmstor.sorted_pts

    lambda = fmmpars.lambda
    ier1 = zeros(Int32,1)
    nlev = tree.nlev
    levelbox = tree.levelbox
    iparentbox = tree.iparentbox
    ichildbox = tree.ichildbox
    icolbox = tree.icolbox
    irowbox = tree.irowbox
    nboxes = tree.nboxes
    nblevel = tree.nblevel
    iboxlev = tree.iboxlev
    istartlev = tree.istartlev
    zll = tree.zll
    blength = tree.blength
    ns = sorted_pts.ns
    srcsort = sorted_pts.srcsort
    isrcladder = sorted_pts.isrcladder

    ifcharge = fmmstor.ifcharge
    chargesort = fmmstor.chargesort
    ifdipole = fmmstor.ifdipole
    dipstrsort = fmmstor.dipstrsort
    dipvecsort = fmmstor.dipvecsort
    ifquad = fmmstor.ifquad
    quadstrsort = fmmstor.quadstrsort
    quadvecsort = fmmstor.quadvecsort
    ifoct = fmmstor.ifoct
    octstrsort = fmmstor.octstrsort
    octvecsort = fmmstor.octvecsort
    isave = fmmstor.isave
    dsave = fmmstor.dsave
    csave = fmmstor.csave

    mtemp, nt = size(targ)
    nt = convert(Int32,nt)

    ifpot1 = boxfmm2d_booltoint32(ifpot)
    ifgrad1 = boxfmm2d_booltoint32(ifgrad)
    ifhess1 = boxfmm2d_booltoint32(ifhess)

    ifder31 = zeros(Int32,1)
    der3targ = zeros(Float64,4,1)

    ccall( (:mbhfmm2d3_targ_,LIBMBHFMM2D), Cvoid,
           (Ref{Float64},Ref{Int32},Ref{Int32},
            Ref{Int32},Ref{Int32},Ref{Int32},Ref{Int32},
            Ref{Int32},Ref{Int32},Ref{Int32},Ref{Int32},
            Ref{Int32},Ref{Float64},
            Ref{Float64},Ref{Int32},Ref{Float64},Ref{Int32},
            Ref{Int32},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Complex{Float64}},
            Ref{Int32},Ref{Float64},Ref{Int32},
            Ref{Float64},Ref{Int32},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Int32},
            Ref{Float64}),
           lambda,ier1,nlev,levelbox,iparentbox,
           ichildbox,icolbox,irowbox,nboxes,nblevel,
           iboxlev,istartlev,zll,
           blength,ns,srcsort,isrcladder,ifcharge,chargesort,
           ifdipole,dipstrsort,dipvecsort,ifquad,
           quadstrsort,quadvecsort,ifoct,octstrsort,
           octvecsort,isave,dsave,csave,nt,targ,ifpot1,
           pottarg,ifgrad1,gradtarg,ifhess1,hesstarg,
           ifder31,der3targ)

    return
end

function mbhfmm2d_targ_exrad!(fmmpars::MBHFMM2DParams,
                              fmmstor::MBHFMM2DStorage,
                              exradsort::Array{Float64,1},
                              targ::Array{Float64,2},
                           ifpot::Bool,pottarg::Array{Float64,1},
                           ifgrad::Bool,gradtarg::Array{Float64,2},
                           ifhess::Bool,hesstarg::Array{Float64,2})

    # unpack

    tree = fmmstor.tree
    sorted_pts = fmmstor.sorted_pts

    lambda = fmmpars.lambda
    ier1 = zeros(Int32,1)
    nlev = tree.nlev
    levelbox = tree.levelbox
    iparentbox = tree.iparentbox
    ichildbox = tree.ichildbox
    icolbox = tree.icolbox
    irowbox = tree.irowbox
    nboxes = tree.nboxes
    nblevel = tree.nblevel
    iboxlev = tree.iboxlev
    istartlev = tree.istartlev
    zll = tree.zll
    blength = tree.blength
    ns = sorted_pts.ns
    srcsort = sorted_pts.srcsort
    isrcladder = sorted_pts.isrcladder

    ifcharge = fmmstor.ifcharge
    chargesort = fmmstor.chargesort
    ifdipole = fmmstor.ifdipole
    dipstrsort = fmmstor.dipstrsort
    dipvecsort = fmmstor.dipvecsort
    ifquad = fmmstor.ifquad
    quadstrsort = fmmstor.quadstrsort
    quadvecsort = fmmstor.quadvecsort
    ifoct = fmmstor.ifoct
    octstrsort = fmmstor.octstrsort
    octvecsort = fmmstor.octvecsort
    isave = fmmstor.isave
    dsave = fmmstor.dsave
    csave = fmmstor.csave

    mtemp, nt = size(targ)
    nt = convert(Int32,nt)

    ifpot1 = boxfmm2d_booltoint32(ifpot)
    ifgrad1 = boxfmm2d_booltoint32(ifgrad)
    ifhess1 = boxfmm2d_booltoint32(ifhess)

    ifder31 = zeros(Int32,1)
    der3targ = zeros(Float64,4,1)

    ccall( (:mbhfmm2d3_targ_exrad_,LIBMBHFMM2D), Cvoid,
           (Ref{Float64},Ref{Int32},Ref{Int32},
            Ref{Int32},Ref{Int32},Ref{Int32},Ref{Int32},
            Ref{Int32},Ref{Int32},Ref{Int32},Ref{Int32},
            Ref{Int32},Ref{Float64},
            Ref{Float64},Ref{Int32},Ref{Float64},Ref{Int32},
            Ref{Float64},Ref{Int32},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Complex{Float64}},
            Ref{Int32},Ref{Float64},Ref{Int32},
            Ref{Float64},Ref{Int32},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Int32},
            Ref{Float64}),
           lambda,ier1,nlev,levelbox,iparentbox,
           ichildbox,icolbox,irowbox,nboxes,nblevel,
           iboxlev,istartlev,zll,
           blength,ns,srcsort,isrcladder,exradsort,
           ifcharge,chargesort,
           ifdipole,dipstrsort,dipvecsort,ifquad,
           quadstrsort,quadvecsort,ifoct,octstrsort,
           octvecsort,isave,dsave,csave,nt,targ,ifpot1,
           pottarg,ifgrad1,gradtarg,ifhess1,hesstarg,
           ifder31,der3targ)

    return
end

function mbhfmm2d_srcsrc!(fmmpars::MBHFMM2DParams,
                        fmmstor::MBHFMM2DStorage,
                        ifpot::Bool,pottarg::Array{Float64,1},
                        ifgrad::Bool,gradtarg::Array{Float64,2},
                        ifhess::Bool,hesstarg::Array{Float64,2})

    # temporary storage (you get the sorted versions from
    # the fmm)

    potsort = copy(pottarg)
    gradsort = copy(gradtarg)
    hesssort = copy(hesstarg)

    # unpack

    tree = fmmstor.tree
    sorted_pts = fmmstor.sorted_pts

    lambda = fmmpars.lambda
    ier1 = zeros(Int32,1)
    nlev = tree.nlev
    levelbox = tree.levelbox
    iparentbox = tree.iparentbox
    ichildbox = tree.ichildbox
    icolbox = tree.icolbox
    irowbox = tree.irowbox
    nboxes = tree.nboxes
    nblevel = tree.nblevel
    iboxlev = tree.iboxlev
    istartlev = tree.istartlev
    zll = tree.zll
    blength = tree.blength
    ns = sorted_pts.ns
    srcsort = sorted_pts.srcsort
    isrcladder = sorted_pts.isrcladder
    isrcsort = sorted_pts.isrcsort

    ifcharge = fmmstor.ifcharge
    chargesort = fmmstor.chargesort
    ifdipole = fmmstor.ifdipole
    dipstrsort = fmmstor.dipstrsort
    dipvecsort = fmmstor.dipvecsort
    ifquad = fmmstor.ifquad
    quadstrsort = fmmstor.quadstrsort
    quadvecsort = fmmstor.quadvecsort
    ifoct = fmmstor.ifoct
    octstrsort = fmmstor.octstrsort
    octvecsort = fmmstor.octvecsort
    isave = fmmstor.isave
    dsave = fmmstor.dsave
    csave = fmmstor.csave

    ifpot1 = boxfmm2d_booltoint32(ifpot)
    ifgrad1 = boxfmm2d_booltoint32(ifgrad)
    ifhess1 = boxfmm2d_booltoint32(ifhess)

    ifder31 = zeros(Int32,1)
    der3targ = zeros(Float64,4,1)

    ccall( (:mbhfmm2d3_srcsrc_,LIBMBHFMM2D), Cvoid,
           (Ref{Float64},Ref{Int32},Ref{Int32},
            Ref{Int32},Ref{Int32},Ref{Int32},Ref{Int32},
            Ref{Int32},Ref{Int32},Ref{Int32},Ref{Int32},
            Ref{Int32},Ref{Float64},
            Ref{Float64},Ref{Int32},Ref{Float64},Ref{Int32},
            Ref{Int32},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Float64},
            Ref{Int32},Ref{Float64},Ref{Complex{Float64}},
            Ref{Int32},Ref{Float64},Ref{Int32},
            Ref{Float64},Ref{Int32},Ref{Float64},
            Ref{Int32},Ref{Float64}),
           lambda,ier1,nlev,levelbox,iparentbox,
           ichildbox,icolbox,irowbox,nboxes,nblevel,
           iboxlev,istartlev,zll,
           blength,ns,srcsort,isrcladder,ifcharge,chargesort,
           ifdipole,dipstrsort,dipvecsort,ifquad,
           quadstrsort,quadvecsort,ifoct,octstrsort,
           octvecsort,isave,dsave,csave,ifpot1,
           potsort,ifgrad1,gradsort,ifhess1,hesssort,
           ifder31,der3targ)

    if ifpot
        for i = 1:ns
            pottarg[isrcsort[i]] = potsort[i]
        end
    end
    if ifgrad
        for i = 1:ns
            gradtarg[:,isrcsort[i]] = gradsort[:,i]
        end
    end
    if ifhess
        for i = 1:ns
            hesstarg[:,isrcsort[i]] = hesssort[:,i]
        end
    end

    return
end

function mbhfmm2d_direct!(fmmpars::MBHFMM2DParams,
                          targ::Array{Float64,2},
                          ifpot::Bool,pottarg::Array{Float64,1},
                          ifgrad::Bool,gradtarg::Array{Float64,2},
                          ifhess::Bool,hesstarg::Array{Float64,2};
                          exrad::Array{Float64,1}=Array{Float64}(undef,0))

    src = fmmpars.src
    mtemp, ns = size(src)
    ns = convert(Int32,ns)
    mtemp, nt = size(targ)
    nt = convert(Int32,nt)

    ifexrad = false
    if length(exrad) > 0
        @assert length(exrad) == ns "length of exrad should equal number of sources"
        ifexrad = true
    end

    lambda = fmmpars.lambda
    ifcharge = boxfmm2d_booltoint32(fmmpars.ifcharge)
    ifdipole = boxfmm2d_booltoint32(fmmpars.ifdipole)
    ifquad = boxfmm2d_booltoint32(fmmpars.ifquad)
    ifoct = boxfmm2d_booltoint32(fmmpars.ifoct)

    charge = fmmpars.charge
    dipstr = fmmpars.dipstr
    dipvec = fmmpars.dipvec
    quadstr = fmmpars.quadstr
    quadvec = fmmpars.quadvec
    octstr = fmmpars.octstr
    octvec = fmmpars.octvec

    ifpot1 = boxfmm2d_booltoint32(ifpot)
    ifgrad1 = boxfmm2d_booltoint32(ifgrad)
    ifhess1 = boxfmm2d_booltoint32(ifhess)

    pottemp = zeros(Float64,1)
    gradtemp = zeros(Float64,2)
    hesstemp = zeros(Float64,3)
    der3temp = zeros(Float64,4)

    if ifexrad
        for i = 1:nt
            pottemp = zeros(Float64,1)
            gradtemp = zeros(Float64,2)
            hesstemp = zeros(Float64,3)
            der3temp = zeros(Float64,4)

            ccall((:mbhpotgrad2dall_cdqo3_add_exrad_,LIBMBHFMM2D),Cvoid,
                  (Ref{Float64},Ref{Float64},Ref{Float64},
                   Ref{Int32},Ref{Int32},
                   Ref{Float64},Ref{Int32},Ref{Float64},
                   Ref{Float64},
                   Ref{Int32},Ref{Float64},Ref{Float64},
                   Ref{Int32},Ref{Float64},Ref{Float64},
                   Ref{Float64},Ref{Int32},Ref{Float64},
                   Ref{Int32},Ref{Float64},Ref{Int32},
                   Ref{Float64}),
                  lambda,src,exrad,ns,ifcharge,charge,
                  ifdipole,dipstr,dipvec,
                  ifquad,quadstr,quadvec,ifoct,octstr,octvec,
                  targ[:,i],ifpot1,pottemp,ifgrad1,gradtemp,
                  ifhess1,hesstemp)

            pottarg[i] = pottemp[1]
            gradtarg[:,i] = gradtemp
            hesstarg[:,i] = hesstemp

        end
    else
        for i = 1:nt

            ccall((:mbhpotgrad2dall_cdqo_,LIBMBHFMM2D),Cvoid,
                  (Ref{Float64},Ref{Float64},Ref{Int32},Ref{Int32},
                   Ref{Float64},Ref{Int32},Ref{Float64},
                   Ref{Float64},
                   Ref{Int32},Ref{Float64},Ref{Float64},
                   Ref{Int32},Ref{Float64},Ref{Float64},
                   Ref{Float64},Ref{Int32},Ref{Float64},
                   Ref{Int32},Ref{Float64},Ref{Int32},
                   Ref{Float64}),
                  lambda,src,ns,ifcharge,charge,ifdipole,
                  dipstr,dipvec,
                  ifquad,quadstr,quadvec,ifoct,octstr,octvec,
                  targ[:,i],ifpot1,pottemp,ifgrad1,gradtemp,
                  ifhess1,hesstemp)

            pottarg[i] = pottemp[1]
            gradtarg[:,i] = gradtemp
            hesstarg[:,i] = hesstemp

        end
    end
    return
end

function mbhfmm2d_direct_self!(fmmpars::MBHFMM2DParams,
                          ifpot::Bool,pottarg::Array{Float64,1},
                          ifgrad::Bool,gradtarg::Array{Float64,2},
                          ifhess::Bool,hesstarg::Array{Float64,2})

    src = fmmpars.src
    mtemp, ns = size(src)
    ns = convert(Int32,ns)

    lambda = fmmpars.lambda
    ifcharge = boxfmm2d_booltoint32(fmmpars.ifcharge)
    ifdipole = boxfmm2d_booltoint32(fmmpars.ifdipole)
    ifquad = boxfmm2d_booltoint32(fmmpars.ifquad)
    ifoct = boxfmm2d_booltoint32(fmmpars.ifoct)

    charge = fmmpars.charge
    dipstr = fmmpars.dipstr
    dipvec = fmmpars.dipvec
    quadstr = fmmpars.quadstr
    quadvec = fmmpars.quadvec
    octstr = fmmpars.octstr
    octvec = fmmpars.octvec

    ifpot1 = boxfmm2d_booltoint32(ifpot)
    ifgrad1 = boxfmm2d_booltoint32(ifgrad)
    ifhess1 = boxfmm2d_booltoint32(ifhess)

    ifder31 = zeros(Int32,1)
    der3targ = zeros(Float64,4,1)

    ccall((:mbhpotgrad2dall_cdqo3_self_,LIBMBHFMM2D),Cvoid,
          (Ref{Float64},Ref{Float64},Ref{Int32},Ref{Int32},
           Ref{Float64},Ref{Int32},Ref{Float64},Ref{Float64},
           Ref{Int32},Ref{Float64},Ref{Float64},
           Ref{Int32},Ref{Float64},Ref{Float64},
           Ref{Int32},Ref{Float64},Ref{Int32},Ref{Float64},
           Ref{Int32},Ref{Float64},Ref{Int32},Ref{Float64}),
          lambda,src,ns,ifcharge,charge,ifdipole,dipstr,dipvec,
          ifquad,quadstr,quadvec,ifoct,octstr,octvec,
          ifpot1,pottarg,ifgrad1,gradtarg,
          ifhess1,hesstarg,ifder31,der3targ)

    return
end