text:
#unum-convert.jl
#implements conversions between unums and ints, floats.

################################################################################
## UNUM TO UNUM conversion

function __subnormal_ubit_trim{ESS,FSS}(::Type{Unum{ESS,FSS}}, esize::UInt16, fsize::UInt16)
  #this function handles the situation where we're trying to convert a zero ulp
  #into a zero ulp of another size, which may or may not have the appropriate
  #capacity to handle the ulp.
  #calculate the (negative) exponent + 1 on the high value of the subnormal number
  high_value_rep::UInt16 = (0x0001 << esize) + fsize
  #find the appropriate esize to rerepresent this value.

  suggested_esize::UInt16 = 0x000F - leading_zeros(high_value_rep)
  esize::UInt16 = min(suggested_esize, max_esize(ESS))

  suggested_fsize::UInt16 = high_value_rep - (0x0001 << esize)
  fsize::UInt16 = min(suggested_fsize, max_fsize(FSS))

  (esize, fsize)
end

function check_lower_cells{ESS,FSS, DEST_FSS}(x::Unum{ESS,FSS}, ::Type{Val{DEST_FSS}})
  if (DEST_FSS < 7)
    start_idx = 2
    accum = (FSS < 7) ? (x.fraction & mask_bot(DEST_FSS)) : (x.fraction.a[1] & mask_bot(DEST_FSS))
  else
    (__cell_length(DEST_FSS) + 1)
    accum = z64
  end

  for idx = start_idx:__cell_length(FSS)
    accum |= x.fraction.a[idx]
  end

  UNUM_UBIT_MASK * (accum != z64)
end

#int the case that DEST_FSS is strictly less than src FSS, then the opration is to pull the single
#value (if DEST_FSS requires an int64)
function pull_upper_cells{ESS,FSS, DEST_FSS}(src::Unum{ESS,FSS}, ::Type{Val{DEST_FSS}})
  if (DEST_FSS < 7)
    return (FSS < 7) ? src.fraction : src.fraction.a[1]
  else
    return ArrayNum{DEST_FSS}(src.fraction.a[1:__cell_length(DEST_FSS)])
  end
end

#in the case that the DEST_FSS is bigger than or equal to the src FSS, then the relevant operation
#is to simply copy over the cells into a blank array with zeros in the rest of the slots.
function copy_cells{ESS, FSS, DEST_FSS}(src::Unum{ESS,FSS}, ::Type{Val{DEST_FSS}})
  (DEST_FSS < 7) && return src.fraction
  result = zero(ArrayNum{DEST_FSS})
  (FSS < 7) ? (result.a[1] = src.fraction) : (result.a[1:__cell_length(FSS)] = src.fraction.a)
  return result
end

doc"""
  `Unums.full_decode(::Unum,::Type{Val{DEST_FSS}})` decodes a unum value into a
  normalized int64 exponent, fraction, fsize, and ubit.  This can be a costly
  operation, since we don't expect too many conversions between unum types.
"""
function full_decode{ESS, FSS, DEST_FSS}(x::Unum{ESS, FSS}, ::Type{Val{DEST_FSS}})
  temp = copy(x)
  leftshift::UInt16 = z16
  subnormal_adjustment::Int64 = zero(Int64)  #this is the extra value we add in for subnormal numbers.

  if (temp.exponent == 0) #then we have a strictly subnormal (strange or regular) number.  Be prepared to shift left.
    #now, count leading zeros.
    leftshift = clz(temp.fraction) + o16
    #next, shift the shadow fraction to the left appropriately.
    frac_lsh!(temp, leftshift)
    subnormal_adjustment = one(Int64)
  end

  #set exponent and fsize, and a default ubit value.
  exponent = decode_exp(temp) - leftshift + subnormal_adjustment

  fsize = min(temp.fsize - leftshift, max_fsize(DEST_FSS))
  ubit::UInt16 = z16

  #reshape the fraction as necessary.
  (DEST_FSS < FSS) && (ubit |= check_lower_cells(temp, Val{DEST_FSS}))

  fraction = (DEST_FSS < FSS) ? pull_upper_cells(temp, Val{DEST_FSS}) : copy_cells(temp, Val{DEST_FSS})

  #full_decode returns the exponent and fraction
  return (exponent, fraction, fsize, ubit)
end

function Base.convert{DEST_ESS,DEST_FSS,SRC_ESS,SRC_FSS}(::Type{Unum{DEST_ESS,DEST_FSS}}, x::Unum{SRC_ESS,SRC_FSS})
  (DEST_ESS == SRC_ESS) && (DEST_FSS == SRC_FSS) && throw(ArgumentError("error attempting to convert the same Unum"))

  #set the type.
  T = (DEST_FSS < 7) ? UnumSmall{DEST_ESS, DEST_FSS} : UnumLarge{DEST_ESS, DEST_FSS}
  A = (DEST_FSS < 7) ? UInt64 : ArrayNum{DEST_FSS}

  #deal with NaN and inf, because those have unusual rules.
  is_nan(x) && return nan(T)
  is_inf(x) && return inf(T, @signof(x))
  #check zeroish numbers (zero ulps and zeros)
  if (x.exponent == z64) && (is_all_zero(x.fraction))
    is_exact(x) && return zero(T)
    #if we're a zero ulp, then we need to do a special conversion done by indexing
    #zero ulps.
    zulp_index = zero_ulp_index(x.esize, x.fsize)
    (esize, fsize) = zero_ulp_params(T, zulp_index)
    return T(z64, zero(A), x.flags, esize, fsize)
  end

  #first convert to a "universal" format that doesn't have restrictions, or
  #subnormals).  be sure to pass the DEST_FSS, so true_fraction can be the
  #needed type

  (true_exponent, true_fraction, true_fsize, ubit) = full_decode(x, Val{DEST_FSS})

  #first, do the exponent part..
  (SRC_ESS <= DEST_ESS) && is_mmr(x) && return nan(T)

  #use the buildunum function to build the unum
  buildunum(T, true_exponent, true_fraction, x.flags | ubit, true_fsize)
end

################################################################################
# convert trampoline: splits a generic unum convert into a specific convert.
function Base.convert{ESS,FSS}(::Type{Unum{ESS,FSS}}, x::Integer)
  (FSS < 7) ? convert(UnumSmall{ESS,FSS}, x) : convert(UnumLarge{ESS,FSS}, x)
end
##################################################################
## INTEGER TO UNUM

const __unsupported_int_types = [BigInt, Int128, UInt128]
#CONVERSIONS - INTEGER -> UNUM
@universal function Base.convert(T::Type{Unum}, i::Integer)
  #currently unsupported:  BigInt, Int128, UInt128
  (typeof(i) in __unsupported_int_types) && throw(ArgumentError("conversion from int type $(typeof(x)) still unsupported"))

  i64val::UInt64 = zero(UInt64)

  #in ESS = 0 we are required to use subnormal one, so this requires
  #special code.
  if (ESS == 0)
    (i == 1) && return one(U)
  end
  #do a zero check
  if (i == 0)
    return zero(U)
  elseif (i < 0)
    #flip the sign and promote the integer to Unt64
    i64val = UInt64(-i)
    flags = UNUM_SIGN_MASK
  else
    #promote to UInt64
    i64val = UInt64(i)
    flags = z16
  end

  #find the msb of x, this will tell us how much to move things
  msbx = 63 - leading_zeros(i64val)
  if (FSS < 7)
    frac = i64val << (64 - msbx)
  else
    frac = zero(ArrayNum{FSS})
    frac[1] = i64val << (64 - msbx)
  end
  fsize = 0x003f  #set it to 63

  r = buildunum(U, msbx, frac, flags, fsize)

  #check for the "infinity hack" where we accidentally generate infinity by having
  #just the right set of bits.
  res = is_inf(r) ? mmr(U, flags & UNUM_SIGN_MASK) : r

  res
end

##################################################################
## FLOATING POINT CONVERSIONS

#create a type for floating point properties
immutable FProp
  I::Type
  ESS::Int
  FSS::Int
  esize::UInt16
  fsize::UInt16
end

#store floating point properties in a dict
const __fp_props = Dict{Type,FProp}(
  Float16 => FProp(UInt16, 3, 4, UInt16(4),  UInt16(9)),
  Float32 => FProp(UInt32, 4, 5, UInt16(7),  UInt16(22)),
  Float64 => FProp(UInt64, 4, 6, UInt16(10), UInt16(51)))

##################################################################
## FLOATS TO UNUM
macro ccreate()
  code = :()
  for F in [Float16, Float32, Float64]
    code = quote
      $code

      function default_convert(x::$F)
        props = __fp_props[$F]
        #generate some bit masks & corresponding shifts
        signbit = (one(UInt64) << (props.esize + props.fsize + 2))
        signshift = (props.esize + props.fsize + 1)
        exponentbits = (signbit - (1 << (props.fsize + 1)))
        exponentshift = props.fsize + 1
        fractionbits = (one(UInt64) << (props.fsize + 1)) - 1
        fractionshift = 64 - props.fsize - 1
        i::UInt64 = reinterpret(props.I, x)
        flags::UInt16 = (i & signbit) >> (signshift)
        exponent = (i & exponentbits) >> (exponentshift)
        fraction = (i & fractionbits) << (fractionshift)
        isnan(x) && return nan(UnumSmall{props.ESS,props.FSS})
        isinf(x) && return inf(UnumSmall{props.ESS,props.FSS}, flags)
        UnumSmall{props.ESS,props.FSS}(exponent, fraction, flags, props.esize, props.fsize)
      end
    end
  end
  esc(code)
end

@ccreate

export default_convert

doc"""
`default_convert` takes floating point numbers and converts them to the equivalent
unums, using the trivial bitshifiting transformation.
"""
default_convert

#helper function to convert from different floating point types, using the default unum as an intermediate
function Base.convert{ESS,FSS}(T::Type{Unum{ESS,FSS}}, x::AbstractFloat)
  prop = __fp_props[typeof(x)]
  (ESS == prop.ESS) && (FSS == prop.FSS) && return default_convert(x)
  convert(Unum{ESS,FSS}, default_convert(x))
end
@universal function Base.convert(T::Type{Unum}, x::AbstractFloat)
  prop = __fp_props[typeof(x)]
  (ESS == prop.ESS) && (FSS == prop.FSS) && return default_convert(x)
  convert(Unum{ESS,FSS}, default_convert(x))
end

##################################################################
## UNUMS TO FLOAT
macro fcreate()
  code = :()
  topbits = Dict(Float16 => 0x8000, Float32 => 0x8000_0000, Float64 => 0x8000_0000_0000_0000)
  for F in [Float16, Float32, Float64]
    topbit = topbits[F]
    code = quote
      $code
      @universal function Base.convert(T::Type{$F}, x::Unum)
        #easy cases
        isnan(x) && return convert(T, NaN)
        is_pos_inf(x) && return convert(T, Inf)
        is_neg_inf(x) && return -convert(T, Inf)
        is_zero(x) && return zero(T)

        fp = __fp_props[$F]
        bits = fp.esize + fp.fsize + 1     #how many total bits
        ebias = 2 ^ (fp.esize) - 1     #exponent bias (= _emax)
        emin = -(ebias) + 1                #minimum exponent
        #create a dummy value that will hold our result.
        res = zero(fp.I)

        signshift = bits - 2

        #highjack the full_decode function to retrieve the unbiased exponent and the fraction.
        (unbiased_exp, src_frac, _, __) = full_decode(x, Val{6})

        #check to see that unbiased_exp is within appropriate bounds for Float32
        (unbiased_exp > ebias) && return convert($F, Inf) * ((x.flags & UNUM_SIGN_MASK == 0) ? 1 : -1)

        if (unbiased_exp < emin)
         delta = emin - unbiased_exp
         unbiased_exp = emin
         fraction = src_frac >> delta
        else
         fraction = src_frac
        end

        #calculate the rebiased exponent and push it into the result.
        res |= convert(fp.I, unbiased_exp + ebias) << (fp.fsize + 1)

        #transfer the fraction bits.
        res |= convert(fp.I, fraction >> (63 - fp.fsize))

        #transfer the signbit
        res |= $topbit * is_negative(x)

        reinterpret(T, res)
      end
    end
  end
  code
end

@fcreate

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
"""
    DAMMfit(Ind_var, Resp, poro_val)

fit the DAMM model parameters to data.

# Examples
```julia-repl
julia> Tₛ = [19.0, 22.0] # soil temperature [°C]
julia> θ = [0.35, 0.22] # soil moisture [m³ m⁻³]
julia> Rₛ = [2.2, 2.8] # respiration observation
julia> Ind_var = hcat(Tₛ, θ)
julia> p = DAMMfit(Ind_var, Rₛ, 0.4)
  10e8
  64
  9.6e-7
  1.3e-4
  0.4
  0.02
  1.0
julia> DAMM(Ind_var, p)
  2.2
  2.8
```
"""
function DAMMfit(Ind_var, Rₛ, poro_val)
  lb = [0.0, 0.0, 0.0, 0.0] # params can't be negative
  p_ini = [1.0, 6.4, 1.0, 1.0] # initial parameters
  p_fact = (1e9, 1e1, 1e-6, 1e-4, 1.0, 1.0, 1.0) # factor of param, better fit
  DAMMfix = (x,p) -> DAMM(x, p_fact.*(p[1], p[2], p[3], p[4], poro_val, 0.02, 1.0))
  fit_bounds = curve_fit(DAMMfix, Ind_var, Rₛ, p_ini, lower=lb)
  p = fit_bounds.param # fitted param result
  fparam = p_fact[1:4].*p # multiply them by p_fact
  params = (fparam[1], fparam[2], fparam[3], fparam[4], poro_val, 0.02, 1.0)
  return params
end

ext: jl


text:
using HullReachNV
using Test
using LazySets

@testset "HullReachNV.jl" begin
    nnet = read_nnet("nnet/test_nnet.nnet")

    solver = HullReach(0.001, false)

    center = [0.5, 0.5]
    radius = [0.5, 0.5]

    in_hyper = Hyperrectangle(center, radius)

    lower = [-3.5, -1.5]
    upper = [-1.5, 2.0]

    out_hyper = Hyperrectangle(low=lower, high=upper)

    problem = Problem(nnet, in_hyper, out_hyper)

    result = solve(solver, problem)

    @test result.status == :violated

end


text:
# Julia wrapper for header: iio.h
# Automatically generated using Clang.jl
const libiio = "libiio"

function iio_create_scan_context(backend, flags)
    ccall((:iio_create_scan_context, libiio), Ptr{iio_scan_context}, (Cstring, UInt32), backend, flags)
end

function iio_scan_context_destroy(ctx)
    ccall((:iio_scan_context_destroy, libiio), Cvoid, (Ptr{iio_scan_context},), ctx)
end

function iio_scan_context_get_info_list(ctx, info)
    ccall((:iio_scan_context_get_info_list, libiio), ssize_t, (Ptr{iio_scan_context}, Ptr{Ptr{Ptr{iio_context_info}}}), ctx, info)
end

function iio_context_info_list_free(info)
    ccall((:iio_context_info_list_free, libiio), Cvoid, (Ptr{Ptr{iio_context_info}},), info)
end

function iio_context_info_get_description(info)
    ccall((:iio_context_info_get_description, libiio), Cstring, (Ptr{iio_context_info},), info)
end

function iio_context_info_get_uri(info)
    ccall((:iio_context_info_get_uri, libiio), Cstring, (Ptr{iio_context_info},), info)
end

function iio_create_scan_block(backend, flags)
    ccall((:iio_create_scan_block, libiio), Ptr{iio_scan_block}, (Cstring, UInt32), backend, flags)
end

function iio_scan_block_destroy(blk)
    ccall((:iio_scan_block_destroy, libiio), Cvoid, (Ptr{iio_scan_block},), blk)
end

function iio_scan_block_scan(blk)
    ccall((:iio_scan_block_scan, libiio), ssize_t, (Ptr{iio_scan_block},), blk)
end

function iio_scan_block_get_info(blk, index)
    ccall((:iio_scan_block_get_info, libiio), Ptr{iio_context_info}, (Ptr{iio_scan_block}, UInt32), blk, index)
end

function iio_library_get_version(major, minor, git_tag)
    ccall((:iio_library_get_version, libiio), Cvoid, (Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt8}), major, minor, git_tag)
end

function iio_strerror(err, dst, len)
    ccall((:iio_strerror, libiio), Cvoid, (Cint, Cstring, Csize_t), err, dst, len)
end

function iio_has_backend(backend)
    ccall((:iio_has_backend, libiio), Bool, (Cstring,), backend)
end

function iio_get_backends_count()
    ccall((:iio_get_backends_count, libiio), UInt32, ())
end

function iio_get_backend(index)
    ccall((:iio_get_backend, libiio), Cstring, (UInt32,), index)
end

function iio_create_default_context()
    ccall((:iio_create_default_context, libiio), Ptr{iio_context}, ())
end

function iio_create_local_context()
    ccall((:iio_create_local_context, libiio), Ptr{iio_context}, ())
end

function iio_create_xml_context(xml_file)
    ccall((:iio_create_xml_context, libiio), Ptr{iio_context}, (Cstring,), xml_file)
end

function iio_create_xml_context_mem(xml, len)
    ccall((:iio_create_xml_context_mem, libiio), Ptr{iio_context}, (Cstring, Csize_t), xml, len)
end

function iio_create_network_context(host)
    ccall((:iio_create_network_context, libiio), Ptr{iio_context}, (Cstring,), host)
end

function iio_create_context_from_uri(uri)
    ccall((:iio_create_context_from_uri, libiio), Ptr{iio_context}, (Cstring,), uri)
end

function iio_context_clone(ctx)
    ccall((:iio_context_clone, libiio), Ptr{iio_context}, (Ptr{iio_context},), ctx)
end

function iio_context_destroy(ctx)
    ccall((:iio_context_destroy, libiio), Cvoid, (Ptr{iio_context},), ctx)
end

function iio_context_get_version(ctx, major, minor, git_tag)
    ccall((:iio_context_get_version, libiio), Cint, (Ptr{iio_context}, Ptr{UInt32}, Ptr{UInt32}, Ptr{UInt8}), ctx, major, minor, git_tag)
end

function iio_context_get_xml(ctx)
    ccall((:iio_context_get_xml, libiio), Cstring, (Ptr{iio_context},), ctx)
end

function iio_context_get_name(ctx)
    ccall((:iio_context_get_name, libiio), Cstring, (Ptr{iio_context},), ctx)
end

function iio_context_get_description(ctx)
    ccall((:iio_context_get_description, libiio), Cstring, (Ptr{iio_context},), ctx)
end

function iio_context_get_attrs_count(ctx)
    ccall((:iio_context_get_attrs_count, libiio), UInt32, (Ptr{iio_context},), ctx)
end

function iio_context_get_attr(ctx, index, name, value)
    ccall((:iio_context_get_attr, libiio), Cint, (Ptr{iio_context}, UInt32, Ptr{Cstring}, Ptr{Cstring}), ctx, index, name, value)
end

function iio_context_get_attr_value(ctx, name)
    ccall((:iio_context_get_attr_value, libiio), Cstring, (Ptr{iio_context}, Cstring), ctx, name)
end

function iio_context_get_devices_count(ctx)
    ccall((:iio_context_get_devices_count, libiio), UInt32, (Ptr{iio_context},), ctx)
end

function iio_context_get_device(ctx, index)
    ccall((:iio_context_get_device, libiio), Ptr{iio_device}, (Ptr{iio_context}, UInt32), ctx, index)
end

function iio_context_find_device(ctx, name)
    ccall((:iio_context_find_device, libiio), Ptr{iio_device}, (Ptr{iio_context}, Cstring), ctx, name)
end

function iio_context_set_timeout(ctx, timeout_ms)
    ccall((:iio_context_set_timeout, libiio), Cint, (Ptr{iio_context}, UInt32), ctx, timeout_ms)
end

function iio_device_get_context(dev)
    ccall((:iio_device_get_context, libiio), Ptr{iio_context}, (Ptr{iio_device},), dev)
end

function iio_device_get_id(dev)
    ccall((:iio_device_get_id, libiio), Cstring, (Ptr{iio_device},), dev)
end

function iio_device_get_name(dev)
    ccall((:iio_device_get_name, libiio), Cstring, (Ptr{iio_device},), dev)
end

function iio_device_get_channels_count(dev)
    ccall((:iio_device_get_channels_count, libiio), UInt32, (Ptr{iio_device},), dev)
end

function iio_device_get_attrs_count(dev)
    ccall((:iio_device_get_attrs_count, libiio), UInt32, (Ptr{iio_device},), dev)
end

function iio_device_get_buffer_attrs_count(dev)
    ccall((:iio_device_get_buffer_attrs_count, libiio), UInt32, (Ptr{iio_device},), dev)
end

function iio_device_get_channel(dev, index)
    ccall((:iio_device_get_channel, libiio), Ptr{iio_channel}, (Ptr{iio_device}, UInt32), dev, index)
end

function iio_device_get_attr(dev, index)
    ccall((:iio_device_get_attr, libiio), Cstring, (Ptr{iio_device}, UInt32), dev, index)
end

function iio_device_get_buffer_attr(dev, index)
    ccall((:iio_device_get_buffer_attr, libiio), Cstring, (Ptr{iio_device}, UInt32), dev, index)
end

function iio_device_find_channel(dev, name, output)
    ccall((:iio_device_find_channel, libiio), Ptr{iio_channel}, (Ptr{iio_device}, Cstring, Bool), dev, name, output)
end

function iio_device_find_attr(dev, name)
    ccall((:iio_device_find_attr, libiio), Cstring, (Ptr{iio_device}, Cstring), dev, name)
end

function iio_device_find_buffer_attr(dev, name)
    ccall((:iio_device_find_buffer_attr, libiio), Cstring, (Ptr{iio_device}, Cstring), dev, name)
end

function iio_device_attr_read(dev, attr, dst, len)
    ccall((:iio_device_attr_read, libiio), ssize_t, (Ptr{iio_device}, Cstring, Cstring, Csize_t), dev, attr, dst, len)
end

function iio_device_attr_read_all(dev, cb, data)
    ccall((:iio_device_attr_read_all, libiio), Cint, (Ptr{iio_device}, Ptr{Cvoid}, Ptr{Cvoid}), dev, cb, data)
end

function iio_device_attr_read_bool(dev, attr, val)
    ccall((:iio_device_attr_read_bool, libiio), Cint, (Ptr{iio_device}, Cstring, Ptr{Bool}), dev, attr, val)
end

function iio_device_attr_read_longlong(dev, attr, val)
    ccall((:iio_device_attr_read_longlong, libiio), Cint, (Ptr{iio_device}, Cstring, Ptr{Clonglong}), dev, attr, val)
end

function iio_device_attr_read_double(dev, attr, val)
    ccall((:iio_device_attr_read_double, libiio), Cint, (Ptr{iio_device}, Cstring, Ptr{Cdouble}), dev, attr, val)
end

function iio_device_attr_write(dev, attr, src)
    ccall((:iio_device_attr_write, libiio), ssize_t, (Ptr{iio_device}, Cstring, Cstring), dev, attr, src)
end

function iio_device_attr_write_raw(dev, attr, src, len)
    ccall((:iio_device_attr_write_raw, libiio), ssize_t, (Ptr{iio_device}, Cstring, Ptr{Cvoid}, Csize_t), dev, attr, src, len)
end

function iio_device_attr_write_all(dev, cb, data)
    ccall((:iio_device_attr_write_all, libiio), Cint, (Ptr{iio_device}, Ptr{Cvoid}, Ptr{Cvoid}), dev, cb, data)
end

function iio_device_attr_write_bool(dev, attr, val)
    ccall((:iio_device_attr_write_bool, libiio), Cint, (Ptr{iio_device}, Cstring, Bool), dev, attr, val)
end

function iio_device_attr_write_longlong(dev, attr, val)
    ccall((:iio_device_attr_write_longlong, libiio), Cint, (Ptr{iio_device}, Cstring, Clonglong), dev, attr, val)
end

function iio_device_attr_write_double(dev, attr, val)
    ccall((:iio_device_attr_write_double, libiio), Cint, (Ptr{iio_device}, Cstring, Cdouble), dev, attr, val)
end

function iio_device_buffer_attr_read(dev, attr, dst, len)
    ccall((:iio_device_buffer_attr_read, libiio), ssize_t, (Ptr{iio_device}, Cstring, Cstring, Csize_t), dev, attr, dst, len)
end

function iio_device_buffer_attr_read_all(dev, cb, data)
    ccall((:iio_device_buffer_attr_read_all, libiio), Cint, (Ptr{iio_device}, Ptr{Cvoid}, Ptr{Cvoid}), dev, cb, data)
end

function iio_device_buffer_attr_read_bool(dev, attr, val)
    ccall((:iio_device_buffer_attr_read_bool, libiio), Cint, (Ptr{iio_device}, Cstring, Ptr{Bool}), dev, attr, val)
end

function iio_device_buffer_attr_read_longlong(dev, attr, val)
    ccall((:iio_device_buffer_attr_read_longlong, libiio), Cint, (Ptr{iio_device}, Cstring, Ptr{Clonglong}), dev, attr, val)
end

function iio_device_buffer_attr_read_double(dev, attr, val)
    ccall((:iio_device_buffer_attr_read_double, libiio), Cint, (Ptr{iio_device}, Cstring, Ptr{Cdouble}), dev, attr, val)
end

function iio_device_buffer_attr_write(dev, attr, src)
    ccall((:iio_device_buffer_attr_write, libiio), ssize_t, (Ptr{iio_device}, Cstring, Cstring), dev, attr, src)
end

function iio_device_buffer_attr_write_raw(dev, attr, src, len)
    ccall((:iio_device_buffer_attr_write_raw, libiio), ssize_t, (Ptr{iio_device}, Cstring, Ptr{Cvoid}, Csize_t), dev, attr, src, len)
end

function iio_device_buffer_attr_write_all(dev, cb, data)
    ccall((:iio_device_buffer_attr_write_all, libiio), Cint, (Ptr{iio_device}, Ptr{Cvoid}, Ptr{Cvoid}), dev, cb, data)
end

function iio_device_buffer_attr_write_bool(dev, attr, val)
    ccall((:iio_device_buffer_attr_write_bool, libiio), Cint, (Ptr{iio_device}, Cstring, Bool), dev, attr, val)
end

function iio_device_buffer_attr_write_longlong(dev, attr, val)
    ccall((:iio_device_buffer_attr_write_longlong, libiio), Cint, (Ptr{iio_device}, Cstring, Clonglong), dev, attr, val)
end

function iio_device_buffer_attr_write_double(dev, attr, val)
    ccall((:iio_device_buffer_attr_write_double, libiio), Cint, (Ptr{iio_device}, Cstring, Cdouble), dev, attr, val)
end

function iio_device_set_data(dev, data)
    ccall((:iio_device_set_data, libiio), Cvoid, (Ptr{iio_device}, Ptr{Cvoid}), dev, data)
end

function iio_device_get_data(dev)
    ccall((:iio_device_get_data, libiio), Ptr{Cvoid}, (Ptr{iio_device},), dev)
end

function iio_device_get_trigger(dev, trigger)
    ccall((:iio_device_get_trigger, libiio), Cint, (Ptr{iio_device}, Ptr{Ptr{iio_device}}), dev, trigger)
end

function iio_device_set_trigger(dev, trigger)
    ccall((:iio_device_set_trigger, libiio), Cint, (Ptr{iio_device}, Ptr{iio_device}), dev, trigger)
end

function iio_device_is_trigger(dev)
    ccall((:iio_device_is_trigger, libiio), Bool, (Ptr{iio_device},), dev)
end

function iio_device_set_kernel_buffers_count(dev, nb_buffers)
    ccall((:iio_device_set_kernel_buffers_count, libiio), Cint, (Ptr{iio_device}, UInt32), dev, nb_buffers)
end

function iio_channel_get_device(chn)
    ccall((:iio_channel_get_device, libiio), Ptr{iio_device}, (Ptr{iio_channel},), chn)
end

function iio_channel_get_id(chn)
    ccall((:iio_channel_get_id, libiio), Cstring, (Ptr{iio_channel},), chn)
end

function iio_channel_get_name(chn)
    ccall((:iio_channel_get_name, libiio), Cstring, (Ptr{iio_channel},), chn)
end

function iio_channel_is_output(chn)
    ccall((:iio_channel_is_output, libiio), Bool, (Ptr{iio_channel},), chn)
end

function iio_channel_is_scan_element(chn)
    ccall((:iio_channel_is_scan_element, libiio), Bool, (Ptr{iio_channel},), chn)
end

function iio_channel_get_attrs_count(chn)
    ccall((:iio_channel_get_attrs_count, libiio), UInt32, (Ptr{iio_channel},), chn)
end

function iio_channel_get_attr(chn, index)
    ccall((:iio_channel_get_attr, libiio), Cstring, (Ptr{iio_channel}, UInt32), chn, index)
end

function iio_channel_find_attr(chn, name)
    ccall((:iio_channel_find_attr, libiio), Cstring, (Ptr{iio_channel}, Cstring), chn, name)
end

function iio_channel_attr_get_filename(chn, attr)
    ccall((:iio_channel_attr_get_filename, libiio), Cstring, (Ptr{iio_channel}, Cstring), chn, attr)
end

function iio_channel_attr_read(chn, attr, dst, len)
    ccall((:iio_channel_attr_read, libiio), ssize_t, (Ptr{iio_channel}, Cstring, Cstring, Csize_t), chn, attr, dst, len)
end

function iio_channel_attr_read_all(chn, cb, data)
    ccall((:iio_channel_attr_read_all, libiio), Cint, (Ptr{iio_channel}, Ptr{Cvoid}, Ptr{Cvoid}), chn, cb, data)
end

function iio_channel_attr_read_bool(chn, attr, val)
    ccall((:iio_channel_attr_read_bool, libiio), Cint, (Ptr{iio_channel}, Cstring, Ptr{Bool}), chn, attr, val)
end

function iio_channel_attr_read_longlong(chn, attr, val)
    ccall((:iio_channel_attr_read_longlong, libiio), Cint, (Ptr{iio_channel}, Cstring, Ptr{Clonglong}), chn, attr, val)
end

function iio_channel_attr_read_double(chn, attr, val)
    ccall((:iio_channel_attr_read_double, libiio), Cint, (Ptr{iio_channel}, Cstring, Ptr{Cdouble}), chn, attr, val)
end

function iio_channel_attr_write(chn, attr, src)
    ccall((:iio_channel_attr_write, libiio), ssize_t, (Ptr{iio_channel}, Cstring, Cstring), chn, attr, src)
end

function iio_channel_attr_write_raw(chn, attr, src, len)
    ccall((:iio_channel_attr_write_raw, libiio), ssize_t, (Ptr{iio_channel}, Cstring, Ptr{Cvoid}, Csize_t), chn, attr, src, len)
end

function iio_channel_attr_write_all(chn, cb, data)
    ccall((:iio_channel_attr_write_all, libiio), Cint, (Ptr{iio_channel}, Ptr{Cvoid}, Ptr{Cvoid}), chn, cb, data)
end

function iio_channel_attr_write_bool(chn, attr, val)
    ccall((:iio_channel_attr_write_bool, libiio), Cint, (Ptr{iio_channel}, Cstring, Bool), chn, attr, val)
end

function iio_channel_attr_write_longlong(chn, attr, val)
    ccall((:iio_channel_attr_write_longlong, libiio), Cint, (Ptr{iio_channel}, Cstring, Clonglong), chn, attr, val)
end

function iio_channel_attr_write_double(chn, attr, val)
    ccall((:iio_channel_attr_write_double, libiio), Cint, (Ptr{iio_channel}, Cstring, Cdouble), chn, attr, val)
end

function iio_channel_enable(chn)
    ccall((:iio_channel_enable, libiio), Cvoid, (Ptr{iio_channel},), chn)
end

function iio_channel_disable(chn)
    ccall((:iio_channel_disable, libiio), Cvoid, (Ptr{iio_channel},), chn)
end

function iio_channel_is_enabled(chn)
    ccall((:iio_channel_is_enabled, libiio), Bool, (Ptr{iio_channel},), chn)
end

function iio_channel_read_raw(chn, buffer, dst, len)
    ccall((:iio_channel_read_raw, libiio), Csize_t, (Ptr{iio_channel}, Ptr{iio_buffer}, Ptr{Cvoid}, Csize_t), chn, buffer, dst, len)
end

function iio_channel_read(chn, buffer, dst, len)
    ccall((:iio_channel_read, libiio), Csize_t, (Ptr{iio_channel}, Ptr{iio_buffer}, Ptr{Cvoid}, Csize_t), chn, buffer, dst, len)
end

function iio_channel_write_raw(chn, buffer, src, len)
    ccall((:iio_channel_write_raw, libiio), Csize_t, (Ptr{iio_channel}, Ptr{iio_buffer}, Ptr{Cvoid}, Csize_t), chn, buffer, src, len)
end

function iio_channel_write(chn, buffer, src, len)
    ccall((:iio_channel_write, libiio), Csize_t, (Ptr{iio_channel}, Ptr{iio_buffer}, Ptr{Cvoid}, Csize_t), chn, buffer, src, len)
end

function iio_channel_set_data(chn, data)
    ccall((:iio_channel_set_data, libiio), Cvoid, (Ptr{iio_channel}, Ptr{Cvoid}), chn, data)
end

function iio_channel_get_data(chn)
    ccall((:iio_channel_get_data, libiio), Ptr{Cvoid}, (Ptr{iio_channel},), chn)
end

function iio_channel_get_type(chn)
    ccall((:iio_channel_get_type, libiio), iio_chan_type, (Ptr{iio_channel},), chn)
end

function iio_channel_get_modifier(chn)
    ccall((:iio_channel_get_modifier, libiio), iio_modifier, (Ptr{iio_channel},), chn)
end

function iio_buffer_get_device(buf)
    ccall((:iio_buffer_get_device, libiio), Ptr{iio_device}, (Ptr{iio_buffer},), buf)
end

function iio_device_create_buffer(dev, samples_count, cyclic)
    ccall((:iio_device_create_buffer, libiio), Ptr{iio_buffer}, (Ptr{iio_device}, Csize_t, Bool), dev, samples_count, cyclic)
end

function iio_buffer_destroy(buf)
    ccall((:iio_buffer_destroy, libiio), Cvoid, (Ptr{iio_buffer},), buf)
end

function iio_buffer_get_poll_fd(buf)
    ccall((:iio_buffer_get_poll_fd, libiio), Cint, (Ptr{iio_buffer},), buf)
end

function iio_buffer_set_blocking_mode(buf, blocking)
    ccall((:iio_buffer_set_blocking_mode, libiio), Cint, (Ptr{iio_buffer}, Bool), buf, blocking)
end

function iio_buffer_refill(buf)
    ccall((:iio_buffer_refill, libiio), ssize_t, (Ptr{iio_buffer},), buf)
end

function iio_buffer_push(buf)
    ccall((:iio_buffer_push, libiio), ssize_t, (Ptr{iio_buffer},), buf)
end

function iio_buffer_push_partial(buf, samples_count)
    ccall((:iio_buffer_push_partial, libiio), ssize_t, (Ptr{iio_buffer}, Csize_t), buf, samples_count)
end

function iio_buffer_cancel(buf)
    ccall((:iio_buffer_cancel, libiio), Cvoid, (Ptr{iio_buffer},), buf)
end

function iio_buffer_start(buf)
    ccall((:iio_buffer_start, libiio), Ptr{Cvoid}, (Ptr{iio_buffer},), buf)
end

function iio_buffer_first(buf, chn)
    ccall((:iio_buffer_first, libiio), Ptr{Cvoid}, (Ptr{iio_buffer}, Ptr{iio_channel}), buf, chn)
end

function iio_buffer_step(buf)
    ccall((:iio_buffer_step, libiio), Cptrdiff_t, (Ptr{iio_buffer},), buf)
end

function iio_buffer_end(buf)
    ccall((:iio_buffer_end, libiio), Ptr{Cvoid}, (Ptr{iio_buffer},), buf)
end

function iio_buffer_foreach_sample(buf, callback, data)
    ccall((:iio_buffer_foreach_sample, libiio), ssize_t, (Ptr{iio_buffer}, Ptr{Cvoid}, Ptr{Cvoid}), buf, callback, data)
end

function iio_buffer_set_data(buf, data)
    ccall((:iio_buffer_set_data, libiio), Cvoid, (Ptr{iio_buffer}, Ptr{Cvoid}), buf, data)
end

function iio_buffer_get_data(buf)
    ccall((:iio_buffer_get_data, libiio), Ptr{Cvoid}, (Ptr{iio_buffer},), buf)
end

function iio_device_get_sample_size(dev)
    ccall((:iio_device_get_sample_size, libiio), ssize_t, (Ptr{iio_device},), dev)
end

function iio_channel_get_index(chn)
    ccall((:iio_channel_get_index, libiio), Clong, (Ptr{iio_channel},), chn)
end

function iio_channel_get_data_format(chn)
    ccall((:iio_channel_get_data_format, libiio), Ptr{iio_data_format}, (Ptr{iio_channel},), chn)
end

function iio_channel_convert(chn, dst, src)
    ccall((:iio_channel_convert, libiio), Cvoid, (Ptr{iio_channel}, Ptr{Cvoid}, Ptr{Cvoid}), chn, dst, src)
end

function iio_channel_convert_inverse(chn, dst, src)
    ccall((:iio_channel_convert_inverse, libiio), Cvoid, (Ptr{iio_channel}, Ptr{Cvoid}, Ptr{Cvoid}), chn, dst, src)
end

function iio_device_get_debug_attrs_count(dev)
    ccall((:iio_device_get_debug_attrs_count, libiio), UInt32, (Ptr{iio_device},), dev)
end

function iio_device_get_debug_attr(dev, index)
    ccall((:iio_device_get_debug_attr, libiio), Cstring, (Ptr{iio_device}, UInt32), dev, index)
end

function iio_device_find_debug_attr(dev, name)
    ccall((:iio_device_find_debug_attr, libiio), Cstring, (Ptr{iio_device}, Cstring), dev, name)
end

function iio_device_debug_attr_read(dev, attr, dst, len)
    ccall((:iio_device_debug_attr_read, libiio), ssize_t, (Ptr{iio_device}, Cstring, Cstring, Csize_t), dev, attr, dst, len)
end

function iio_device_debug_attr_read_all(dev, cb, data)
    ccall((:iio_device_debug_attr_read_all, libiio), Cint, (Ptr{iio_device}, Ptr{Cvoid}, Ptr{Cvoid}), dev, cb, data)
end

function iio_device_debug_attr_write(dev, attr, src)
    ccall((:iio_device_debug_attr_write, libiio), ssize_t, (Ptr{iio_device}, Cstring, Cstring), dev, attr, src)
end

function iio_device_debug_attr_write_raw(dev, attr, src, len)
    ccall((:iio_device_debug_attr_write_raw, libiio), ssize_t, (Ptr{iio_device}, Cstring, Ptr{Cvoid}, Csize_t), dev, attr, src, len)
end

function iio_device_debug_attr_write_all(dev, cb, data)
    ccall((:iio_device_debug_attr_write_all, libiio), Cint, (Ptr{iio_device}, Ptr{Cvoid}, Ptr{Cvoid}), dev, cb, data)
end

function iio_device_debug_attr_read_bool(dev, attr, val)
    ccall((:iio_device_debug_attr_read_bool, libiio), Cint, (Ptr{iio_device}, Cstring, Ptr{Bool}), dev, attr, val)
end

function iio_device_debug_attr_read_longlong(dev, attr, val)
    ccall((:iio_device_debug_attr_read_longlong, libiio), Cint, (Ptr{iio_device}, Cstring, Ptr{Clonglong}), dev, attr, val)
end

function iio_device_debug_attr_read_double(dev, attr, val)
    ccall((:iio_device_debug_attr_read_double, libiio), Cint, (Ptr{iio_device}, Cstring, Ptr{Cdouble}), dev, attr, val)
end

function iio_device_debug_attr_write_bool(dev, attr, val)
    ccall((:iio_device_debug_attr_write_bool, libiio), Cint, (Ptr{iio_device}, Cstring, Bool), dev, attr, val)
end

function iio_device_debug_attr_write_longlong(dev, attr, val)
    ccall((:iio_device_debug_attr_write_longlong, libiio), Cint, (Ptr{iio_device}, Cstring, Clonglong), dev, attr, val)
end

function iio_device_debug_attr_write_double(dev, attr, val)
    ccall((:iio_device_debug_attr_write_double, libiio), Cint, (Ptr{iio_device}, Cstring, Cdouble), dev, attr, val)
end

function iio_device_identify_filename(dev, filename, chn, attr)
    ccall((:iio_device_identify_filename, libiio), Cint, (Ptr{iio_device}, Cstring, Ptr{Ptr{iio_channel}}, Ptr{Cstring}), dev, filename, chn, attr)
end

function iio_device_reg_write(dev, address, value)
    ccall((:iio_device_reg_write, libiio), Cint, (Ptr{iio_device}, UInt32, UInt32), dev, address, value)
end

function iio_device_reg_read(dev, address, value)
    ccall((:iio_device_reg_read, libiio), Cint, (Ptr{iio_device}, UInt32, Ptr{UInt32}), dev, address, value)
end


text:
active_events(events::NamedTuple, events_active::NamedTuple, events_save::Tuple) =
  active_events(events, events_active, NamedTuple{keys(events)}(fill(events_save,length(events))))

active_events(events::NamedTuple, events_active::NamedTuple, events_save::Vector{Pair{Symbol, Tuple{Bool, Bool}}}) =
  active_events(events, events_active, NamedTuple(events_save))

function active_events(events::NamedTuple, events_active::NamedTuple, events_save::NamedTuple)
  ev_names = keys((events))
  for k in (keys(events_active)..., keys(events_save)...)
    @assert k in ev_names "Event $k not found."
  end
  _events_save = merge(NamedTuple{ev_names}(fill((true,true), length(events))), events_save)

  return Tuple(add_event(events[ev], _events_save[ev], ev) for ev in keys(events) if events_active[ev])
end

function add_event(evt::TimeEvent, events_save::Tuple{Bool, Bool}=(true,true), evt_name=nothing)
  tstops = Float64[]
  #scn_func(u, t, integrator) = t in tstops

  function init_time_event(cb,u,t,integrator)
      append!(tstops, evt.condition_func(integrator.sol.prob.p.constants, integrator.sol.prob.tspan))
      tf = integrator.sol.prob.tspan[2]
      [add_tstop!(integrator, tstop) for tstop in tstops if tstop <= tf]
      #[add_tstop!(integrator, tstop) for tstop in tstops]
      cb.condition(u,t,integrator) ? cb.affect!(integrator) : nothing
  end

  DiscreteCallback(
        (u,t,integrator) -> t in tstops,
        (integrator) -> evt_func_wrapper(integrator, evt.affect_func, events_save, evt_name),
        initialize = init_time_event,
        save_positions=(false,false)
  )
end

function add_event(evt::CEvent, events_save::Tuple{Bool, Bool}=(true,true), evt_name=nothing)
  ContinuousCallback(
      evt.condition_func,
      (integrator) -> evt_func_wrapper(integrator, evt.affect_func, events_save, evt_name),
      (integrator) -> nothing,
      save_positions=(false,false)
  )
end

function add_event(evt::StopEvent, events_save::Tuple{Bool, Bool}=(true,false), evt_name=nothing)
  DiscreteCallback(
    evt.condition_func,
    (integrator) -> evt_func_wrapper(integrator, terminate!, events_save, evt_name),
    save_positions=(false,false)
  )
end

function evt_func_wrapper(integrator, evt_func, events_save, evt_name)
  # check saveat values before applying a callback
  save_after_step!(integrator)
  # save timepoint before and after applying a callback
  first(events_save) && save_timepoint!(integrator, :ode_) #affect_func!(integrator, true)
  evt_func(integrator)
  last(events_save) && save_timepoint!(integrator, evt_name)
  reset_dt!(integrator)
end

function reset_dt!(integrator::Sundials.AbstractSundialsIntegrator)
# not implemented
end

function reset_dt!(integrator::DiffEqBase.AbstractODEIntegrator)
  if integrator.t != integrator.sol.prob.tspan[1]  # exclude events at zero
    auto_dt_reset!(integrator)
    set_proposed_dt!(integrator, integrator.dt)
  end
end


text:
module TestMakeParameters
  using BioEnergeticFoodWebs
  using Test

  # Test the keyword interface
  correct_network = [0 1 0 0; 0 0 1 1; 0 0 0 0; 0 0 0 0]
  parameters = model_parameters(correct_network, Z=2.0)
  @test parameters[:Z] == 2.0

  # Test that producers, etc, are identified
  @test parameters[:is_producer][4] == true
  @test parameters[:is_producer][3] == true
  @test parameters[:is_producer][2] == false
  @test parameters[:is_producer][1] == false

  @test parameters[:is_herbivore][4] == false
  @test parameters[:is_herbivore][3] == false
  @test parameters[:is_herbivore][2] == true
  @test parameters[:is_herbivore][1] == false

  # Test the direct interface
  correct_network = [0 1 0 0; 0 0 1 1; 0 0 0 0; 0 0 0 0]
  parameters = BioEnergeticFoodWebs.model_parameters(correct_network, Z=2.0)
  @test parameters[:Z] == 2.0

  # Test that there is an exception if the vertebrates is of the wrong size
  wrong_vert = vec([true true true false false true true true false false])
  @test_throws ErrorException model_parameters(correct_network, vertebrates=wrong_vert)

  # Test that there the vertebrates can be passed
  right_vert = vec([true false true false])
  parameters = model_parameters(correct_network, vertebrates=right_vert)
  @test right_vert == parameters[:vertebrates]

  # Test that there is an exception if the body masses is of the wrong size
  wrong_bs = rand(100)
  @test_throws ErrorException model_parameters(correct_network, bodymass=wrong_bs)

  # Test that there the bodymasses can be passed
  right_bs = rand(4)
  parameters = model_parameters(correct_network, bodymass=right_bs)
  @test right_bs == parameters[:bodymass]

  # Test that the metabolic rates are calculated from bodymass
  p = model_parameters(correct_network)
  p_b = model_parameters(correct_network, bodymass = p[:bodymass])
  @test p[:x] == p_b[:x]

  # Test that there is an exception if the wrong productivity is used
  wrong_pr = :global
  @test_throws ErrorException model_parameters(correct_network, productivity=wrong_pr)

  # Test that there the productivity can be passed
  right_pr = :system
  parameters = model_parameters(correct_network, productivity=right_pr)
  @test right_pr == parameters[:productivity]

  # Test that there is an exception if the wrong rewiring method is used
  wrong_rm = :Rewire
  @test_throws ErrorException model_parameters(correct_network, rewire_method=wrong_rm)

  # Test that the rewire method can be passed
  right_rm = :ADBM
  parameters = model_parameters(correct_network, rewire_method=right_rm)
  @test right_rm == parameters[:rewire_method]

  # Test the adbm_parameters function:
  # Test that the parameters can be passed
  parameters = model_parameters(correct_network, rewire_method=:none)
  test_e = 2.0
  test_a_adbm = 0.02
  test_ai = -0.5
  test_aj = -0.5
  test_b = 0.4
  test_h_adbm = 2.0
  test_hi = 2.0
  test_hj = 2.0
  test_n = 2.0
  test_ni= -0.8
  test_Hmethod = :power
  test_Nmethod = :biomass
  BioEnergeticFoodWebs.adbm_parameters(parameters, test_e, test_a_adbm, test_ai, test_aj, test_b, test_h_adbm, test_hi, test_hj, test_n, test_ni, test_Hmethod, test_Nmethod)
  @test test_e == parameters[:e]
  @test test_a_adbm == parameters[:a_adbm]
  @test test_ai == parameters[:ai]
  @test test_aj == parameters[:aj]
  @test test_b == parameters[:b]
  @test test_h_adbm == parameters[:h_adbm]
  @test test_hi == parameters[:hi]
  @test test_hj == parameters[:hj]
  @test test_n == parameters[:n]
  @test test_ni == parameters[:ni]
  @test test_Hmethod == parameters[:Hmethod]
  @test test_Nmethod == parameters[:Nmethod]

  # Test that there is an exception if the wrong Nmethod is passed
  wrong_Nmethod = :logbiomass
  @test_throws ErrorException model_parameters(correct_network, rewire_method=:ADBM, Nmethod=wrong_Nmethod)

  # Test that there is an exception if the wrong Hmethod is passed
  wrong_Hmethod = :sum
  @test_throws ErrorException model_parameters(correct_network, rewire_method=:ADBM, Hmethod=wrong_Hmethod)

  # Test that a cost matrix (each cost = 1.0) of size (S, S) is added to p
  S = size(correct_network, 1)
  right_cm = ones(S,S)
  @test right_cm == parameters[:costMat]

  # Tests for the gilljam_parameters function:
    # Tests for the internal preference_parameters function:

      # parameters can be passed
  test_cost = .2
  test_preferenceMethod = :specialist
  test_specialistPrefMag = .8
  pref_par = BioEnergeticFoodWebs.preference_parameters(test_cost, test_specialistPrefMag, correct_network, test_preferenceMethod)
  @test pref_par[:cost] == test_cost
  @test pref_par[:preferenceMethod] == test_preferenceMethod
  @test pref_par[:specialistPrefMag] == test_specialistPrefMag

      # jaccard similarity matrix
  network_1 = correct_network #all species have different resources
  pref_par_1 = BioEnergeticFoodWebs.preference_parameters(test_cost, test_specialistPrefMag, network_1, test_preferenceMethod)
  n1_simIndex = [[1,2,3,4],[1,2,3,4],[1,2,3,4],[1,2,3,4]]
  @test n1_simIndex == pref_par_1[:similarity]

  network_2 = [0 0 1 1; 0 0 1 1; 0 0 0 0; 0 0 0 0] #species 1 and 2 have the same diet
  pref_par_2 = BioEnergeticFoodWebs.preference_parameters(test_cost, test_specialistPrefMag, network_2, test_preferenceMethod)
  n2_simIndex = [[1,3,4,2],[2,3,4,1],[1,2,3,4],[1,2,3,4]]
  @test n2_simIndex == pref_par_2[:similarity]

      # cost matrix (each cost = 1.0) of size (S, S) is added to p
  S = size(correct_network, 1)
  right_cm = ones(S,S)
  @test right_cm == pref_par[:costMat]

    # Test the internal get_specialist_preferences function:

      # Test that the preferences returned are empty if test_preferenceMethod = :generalist
  pm = :generalist
  pr = BioEnergeticFoodWebs.preference_parameters(test_cost, test_specialistPrefMag, correct_network, pm)
  right_pref = zeros(Int64, S)
  pref = BioEnergeticFoodWebs.get_specialist_preferences(pr, correct_network)
  @test right_pref == pref

      # Test that the preferences returned are correct when test_preferenceMethod = :specialist
  pm = :specialist
  pr = BioEnergeticFoodWebs.preference_parameters(test_cost, test_specialistPrefMag, correct_network, pm)
  pref = BioEnergeticFoodWebs.get_specialist_preferences(pr, correct_network)
  right_pref_1 = [2,4,0,0]
  right_pref_2 = [2,3,0,0]
  @test (pref == right_pref_1) | (pref == right_pref_2)

    # Test the gilljam_parameters function:

      # exception if the wrong preferenceMethod is passed
  wrong_pm = :gen
  @test_throws ErrorException model_parameters(correct_network, rewire_method=:DS, preferenceMethod=wrong_pm)
  parameters = model_parameters(correct_network, rewire_method = :DS, preferenceMethod = :specialist)
  @test parameters[:similarity] == n1_simIndex
  @test parameters[:specialistPrefMag] == .9
  @test parameters[:extinctions] == Int[]
  @test parameters[:preferenceMethod] == :specialist
  @test parameters[:cost] == .0
  @test parameters[:costMat] == right_cm
  @test (parameters[:specialistPref] == right_pref_1) | (parameters[:specialistPref] == right_pref_2)

  # Test that an empty vector parameters[:extinctions] is returned when using Staniczenko's rewiring method
  extinctions = Int[]
  p = model_parameters(correct_network, rewire_method = :DO)
  @test extinctions == parameters[:extinctions]

end

module TestUpdateParameters
  using BioEnergeticFoodWebs
  using Test

  A = [0 0 1 ; 0 0 0 ; 0 0 0] #species 1 (herbivore) feed on 3 (producer)
  b = [0.2, 0.3, 0.0] #extinction of species 3

  # Test with ADBM rewiring method
  RWmethod = :ADBM
  parameters = model_parameters(A, rewire_method = RWmethod)
  old_p = copy(parameters)
  BioEnergeticFoodWebs.update_rewiring_parameters(parameters, b, 1)

  #check that all links from and to 3 are gone
  @test parameters[:A] == Int.(zero(A))
  #check that the parameters have been updated
  #species 1 was an herbivore ...
  @test BioEnergeticFoodWebs.get_herbivores(old_p) == [true, false, false]
  #...but not anymore
  @test BioEnergeticFoodWebs.get_herbivores(parameters) == [false, false, false]
  # it is still a consumer though
  @test parameters[:is_producer] == old_p[:is_producer] == [false, true, true]
  #preferences have been updated
  @test BioEnergeticFoodWebs.getW_preference(old_p) == float.([0 0 1 ; 0 0 0 ; 0 0 0])
  @test BioEnergeticFoodWebs.getW_preference(parameters) == zero(A)
  #efficiency have been updated
  @test BioEnergeticFoodWebs.get_efficiency(old_p) == float.([0 0 old_p[:e_herbivore] ; 0 0 0 ; 0 0 0])
  @test BioEnergeticFoodWebs.get_efficiency(parameters) == zero(A)

  # Test with Gilljam rewiring method
  A = [0 0 1 0 ; 0 0 1 1 ; 0 0 0 0 ; 0 0 0 0]
  b = [0.2, 0.3, 0.0, 0.5] #extinction of species 3
  RWmethod = :DS
  parameters = model_parameters(A, rewire_method = RWmethod)
  old_p = copy(parameters)
  BioEnergeticFoodWebs.update_rewiring_parameters(parameters, b, 1)

  #check that all links from and to 3 are gone
  @test parameters[:A][:,3] ==  parameters[:A][3,:] == zeros(4)
  @test parameters[:extinctions] == [3]
  @test BioEnergeticFoodWebs.get_herbivores(parameters) == [true, true, false, false]
  @test parameters[:is_producer] == old_p[:is_producer] == [false, false, true, true]
  eff = float.(zero(parameters[:A]))
  eff[parameters[:A] .> 0] .= parameters[:e_herbivore]
  @test BioEnergeticFoodWebs.get_efficiency(parameters) == eff
  pref = float.(parameters[:A])
  @test BioEnergeticFoodWebs.getW_preference(parameters) == pref

  #with specialists
  p_specialists = model_parameters(A, rewire_method = RWmethod, preferenceMethod = :specialist)
  p_specialists[:A] = parameters[:A]
  p_specialists[:extinctions] = parameters[:extinctions]
  BioEnergeticFoodWebs.update_specialist_preference(p_specialists)
  @test p_specialists[:specialistPref] == [4, 4, 0, 0]

  # Test with Staniczenko rewiring method
  A = [0 0 1 0 ; 0 0 1 1 ; 0 0 0 0 ; 0 0 0 0]
  b = [0.2, 0.3, 0.0, 0.5]
  RWmethod = :DO
  parameters = model_parameters(A, rewire_method = RWmethod)
  old_p = copy(parameters)
  BioEnergeticFoodWebs.update_rewiring_parameters(parameters, b, 1)
  #no released prey => no new link
  @test parameters[:A][:,3] ==  parameters[:A][3,:] == zeros(4)
  # 1 is not an herbivore anymore (no resource)
  @test BioEnergeticFoodWebs.get_herbivores(parameters) == [false, true, false, false]
  pref = float.(parameters[:A])
  @test BioEnergeticFoodWebs.getW_preference(parameters) == pref
  eff = float.(zero(parameters[:A]))
  eff[parameters[:A] .> 0] .= parameters[:e_herbivore]
  @test BioEnergeticFoodWebs.get_efficiency(parameters) == eff

end

# Nutrient Intake model - parameters

module TestNIParameters
  using BioEnergeticFoodWebs
  using Test

  # Test that the parameters can be passed
  A = [0 1 1 ; 0 0 0 ; 0 0 0]
  k1 = [0, 0.1, 0.2]
  k2 = [0, 0.2, 0.1]
  turnover_rate = 0.3
  s = [3.0, 4.0]
  content = [0.8, 0.6]
  parameters = model_parameters(A, productivity = :nutrients, K1 = k1, K2 = k2, D = turnover_rate, supply = s, υ = content)
  @test parameters[:K1] == k1
  @test parameters[:K2] == k2
  @test parameters[:D] == turnover_rate
  @test parameters[:supply] == s
  @test parameters[:υ] == content

  # Test that the NP parameters are not returned when the function is called with productivity != :nutrient
  parameters = model_parameters(A, K1 = k1, K2 = k2, D = turnover_rate, supply = s, υ = content)
  @test_throws KeyError parameters[:K1]
  @test_throws KeyError parameters[:K2]
  @test_throws KeyError parameters[:D]
  @test_throws KeyError parameters[:supply]
  @test_throws KeyError parameters[:υ]

  # Test that there is an error when the vectors length is incorrect
  k1 = [0.1, 0.2]
  @test_throws ErrorException model_parameters(A, productivity = :nutrients, K1 = k1)
  k2 = [0.2, 0.1]
  @test_throws ErrorException model_parameters(A, productivity = :nutrients, K2 = k2)
  s = [3.0, 2.0, 4.0]
  @test_throws ErrorException model_parameters(A, productivity = :nutrients, supply = s)
  content = [1.0]
  @test_throws ErrorException model_parameters(A, productivity = :nutrients, υ = content)

  # Producer metabolic rates
  parameters = model_parameters(A, productivity = :nutrients)
  @test parameters[:x] == [0.3141, 0.138, 0.138]

end

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
# -----------------------------------------------------------------------------
#
#
#
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
import Pkg;
Pkg.add(name="Vcov", version="0.4.2", preserve=Pkg.PRESERVE_ALL); # this new version has the fix
Pkg.add(name="FixedEffectModels", version="1.4.2", preserve=Pkg.PRESERVE_ALL); # this new version has the fix
Pkg.add.(["FileIO", "DataFrames", "PairsMacros", "Dates"])

using FileIO, DataFrames, PairsMacros, Dates
using FixedEffectModels
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
# --- load the data and clean it up to adapt it for julia
df_reg = load("data/iv_interact_fe_test.dta") |> DataFrame;
transform!(df_reg, [:group, :Z_range, :group_x_Z] .=> ByRow(x->Int(x)) .=> [:group, :Z_range, :group_x_Z]);
for i in 1:4, j in 1:3
  u = Symbol("group_" * string(i) * "_" * string(j));
  transform!(df_reg, @cols($u = Int.($u) ));
  categorical!(df_reg, u);
end
categorical!(df_reg, :group);
categorical!(df_reg, :Z_range);
categorical!(df_reg, :group_x_Z);
# -----------------------------------------------------------------------------

# -----------------------------------------------------------------------------
# julia code works but F-stat is not inline with old version of reghdfe
try
   eval(Meta.parse("reg(df_reg, @formula(Y ~ fe(unit_id) + (X ~ group&Z_range)), Vcov.robust())"))
catch e
   @error "Something went wrong" exception=(e, catch_backtrace())
end

try
   eval(Meta.parse("reg(df_reg, @formula(Y ~ fe(unit_id) + (X ~ group_x_Z)), Vcov.robust())"))
catch e
   @error "Something went wrong" exception=(e, catch_backtrace())
end
# -----------------------------------------------------------------------------

max_stars_repo_name: eloualiche/snippets

max_issues_repo_name: eloualiche/snippets

max_forks_repo_name: eloualiche/snippets


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

name = "ViennaRNA"
version = v"2.4.18"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/ViennaRNA/ViennaRNA/releases/download/v2.4.18/ViennaRNA-2.4.18.tar.gz",
                  "b276cfed7c3bea4821c8272750b5b22fa6fed614d71a98cfb25eb0095369657c")
]

# Bash recipe for building across all platforms
# Notes:
# - linking problems when using (built-in or JLL pkg) SVM (libsvm) on non-linux targets
# - compilation problems when using GSL on non-linux targets
script = raw"""
cd $WORKSPACE/srcdir
cd ViennaRNA-2.4.18/

# configure script needs bash
sed -i -e '1s|#! /bin/sh|#! /bin/bash|' configure

# configure script fails (only in sandbox) without
# setting ac_cv_func_malloc_0_nonnull=yes ac_cv_func_realloc_0_nonnull=yes
ac_cv_func_malloc_0_nonnull=yes ac_cv_func_realloc_0_nonnull=yes \
    ./configure \
    --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --with-pic --disable-c11 \
    --without-perl --without-python --without-python3 \
    --without-doc --without-tutorial --without-tutorial-pdf \
    --without-gsl --without-svm \
    --without-forester

make -j${nproc}
make install

# create and install a shared library libRNA
#cc -shared -o "${libdir}/libRNA.${dlext}" -Wl,$(flagon --whole-archive) ./src/ViennaRNA/libRNA.a -Wl,$(flagon --no-whole-archive) -lgsl -lgslcblas -lm -lsvm -lgomp -lgsl -lmpfr
cc -shared -o "${libdir}/libRNA.${dlext}" -Wl,$(flagon --whole-archive) ./src/ViennaRNA/libRNA.a -Wl,$(flagon --no-whole-archive) -lm -lgomp -lmpfr
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
# remedy for std::string incompatibilities across the GCC 4/5 version boundary
# (suggested by BinaryBuilder)
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("RNAPKplex", :RNAPKplex),
    ExecutableProduct("RNAlocmin", :RNAlocmin),
    ExecutableProduct("RNAplex", :RNAplex),
    ExecutableProduct("RNAeval", :RNAeval),
    ExecutableProduct("RNAsnoop", :RNAsnoop),
    ExecutableProduct("ct2db", :ct2db),
    ExecutableProduct("RNAcofold", :RNAcofold),
    ExecutableProduct("RNApaln", :RNApaln),
    LibraryProduct("libRNA", :libRNA),
    ExecutableProduct("Kinfold", :Kinfold),
    ExecutableProduct("RNApdist", :RNApdist),
    ExecutableProduct("RNAplot", :RNAplot),
    ExecutableProduct("RNApvmin", :RNApvmin),
    ExecutableProduct("RNAaliduplex", :RNAaliduplex),
    ExecutableProduct("RNAduplex", :RNAduplex),
    ExecutableProduct("RNAdistance", :RNAdistance),
    ExecutableProduct("RNAheat", :RNAheat),
    ExecutableProduct("RNALalifold", :RNALalifold),
    ExecutableProduct("popt", :popt),
    ExecutableProduct("RNAalifold", :RNAalifold),
    ExecutableProduct("RNAsubopt", :RNAsubopt),
    ExecutableProduct("RNALfold", :RNALfold),
    ExecutableProduct("RNAfold", :RNAfold),
    ExecutableProduct("RNAinverse", :RNAinverse),
    ExecutableProduct("RNA2Dfold", :RNA2Dfold),
    ExecutableProduct("RNAdos", :RNAdos),
    ExecutableProduct("RNAparconv", :RNAparconv),
    ExecutableProduct("RNAplfold", :RNAplfold),
    ExecutableProduct("RNAup", :RNAup),
    ExecutableProduct("b2ct", :b2ct)
]

# Dependencies that must be installed before this package can be built
# CompilerSupportLibraries_jll suggested by BinaryBuilder for OpenMP support
dependencies = [
    Dependency(PackageSpec(name="MPFR_jll", uuid="3a97d323-0669-5f0c-9066-3539efd106a3"))
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
#    Dependency(PackageSpec(name="GSL_jll", uuid="1b77fbbe-d8ee-58f0-85f9-836ddc23a7a4"))
#    Dependency(PackageSpec(name="libsvm_jll", uuid="08558c22-525a-5d2a-acf6-0ac6658ffce4"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat="1.6", preferred_gcc_version = v"11")

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
@testset "tagmap" begin
    # Visit every entry in the map.
    @testset "tag $t" for t in keys(PNML.tagmap)
        @test haskey(PNML.tagmap,t)
        @test !isempty(methods(PNML.tagmap[t], (EzXML.Node,))) ||
              !isempty(methods(PNML.tagmap[t], (EzXML.Node, PnmlType)))
    end
    #TODO: Add non-trivial tests.
end

@testset "pnmltype" begin
   #@test_call PNML.PnmlTypes.default_pntd_map()
   @test_call PNML.pnmltype(PnmlCore())
   @test_call pntd_symbol("foo")
   @test_call PNML.pnmltype("pnmlcore")
   @test_call PNML.pnmltype(:pnmlcore)
end


text:
module Logging2

using Logging

include("LineBufferedIO.jl")
include("LoggingStream.jl")

#-------------------------------------------------------------------------------
# Utilities

for (redirect_func, stream_name) in [
        (:redirect_stdout, :stdout),
        (:redirect_stderr, :stderr)
    ]
    @eval function Base.$redirect_func(f::Function, logger::AbstractLogger; level=Logging.Info)
        result = nothing
        prev_stream = $stream_name
        output = LineBufferedIO(LoggingStream(logger, id=$(QuoteNode(stream_name)); level=level))
        rd,rw = $redirect_func()
        try
            @sync begin
                try
                    Threads.@spawn write(output, rd) # loops until !eof(rd)
                    result = f()
                finally
                    # To close the read side of the pipe, we must close *all*
                    # writers. This includes `rw`, but *also* the dup'd fd
                    # created behind the scenes by redirect_func(). (To close
                    # that, must call redirect_func() here with the prev stream.)
                    close(rw)
                    $redirect_func(prev_stream)
                end
            end
        finally
            close(rd)
            close(output)
        end
        return result
    end
end

@doc """
    redirect_stdout(f::Function, logger::AbstractLogger)

Redirect the global stdout stream to `logger`, with each line becoming a log
event during the execution of the function `f`.

!!! note
    In contrast to the dynamic scope of the usual logging system frontend (`@info`,
    etc), `stdout` is a global object so it's not entirely clear that we can
    collect the logger from the current dynamic scope where `Base.stdout` is
    looked up, and efficiently use it.

    In particular, some particular uses of stdout require it to have an
    operating system primitive like a `Pipe` backing the object. However not
    all uses require this, and it may be possible to improve the situation in
    the future.

# Examples

Here's how you use `redirect_stdout` in structured concurrency style:

```
redirect_stdout(current_logger()) do
    println("Hi")
    run(`ls`)
end
```
""" redirect_stdout

@doc """
    redirect_stderr(f::Function, logger::AbstractLogger)

Redirect the global stderr stream to `logger`, with each line becoming a log
event during the execution of the function `f`.

See [`redirect_stdout`](@ref) for examples and additional information.
""" redirect_stderr

end


text:
# Schrodinger.jl Tests
using Schrodinger
using Compat, Compat.Test, Compat.LinearAlgebra, Compat.SparseArrays

@testset "All Tests" begin
include("quobj.jl")
include("math.jl")
include("special.jl")
include("dynamics.jl")
include("tomography.jl")
include("misc.jl")
end


text:
isdefined(:LinearInterpolators) || include("../src/LinearInterpolators.jl")

module LinearInterpolatorsKernelsDemos

using LinearInterpolators.Kernels

using PyCall
# pygui(:gtk); # can be: :wx, :gtk, :qt
using LaTeXStrings
import PyPlot
const plt = PyPlot

box = RectangularSpline()
triangle = LinearSpline()
quadratic = QuadraticSpline()
cubic = CubicSpline()
catmull_rom = CatmullRomSpline()
mitchell_netravali = MitchellNetravaliSpline()
keys = KeysSpline(0.5)
lanczos2 = LanczosKernel(2)
lanczos4 = LanczosKernel(4)
lanczos6 = LanczosKernel(6)
lanczos8 = LanczosKernel(8)
lanczos10 = LanczosKernel(10)

plt.figure(10)
plt.clf()
x = linspace(-6,7,1000)
plt.plot(x, box(x), color="darkgreen",
         linewidth=2.0, linestyle="-", label="box")
plt.plot(x, triangle(x), color="darkblue",
         linewidth=2.0, linestyle="-", label="triangle (linear B-spline)")
plt.plot(x, quadratic(x), color="darkcyan",
         linewidth=2.0, linestyle="-", label="quadratic B-spline")
plt.plot(x, cubic(x), color="firebrick",
         linewidth=2.0, linestyle="-", label="cubic B-spline")
plt.plot(x, catmull_rom(x), color="orange",
         linewidth=2.0, linestyle="-", label="Catmull-Rom spline")
plt.plot(x, mitchell_netravali(x), color="violet",
         linewidth=2.0, linestyle="-", label="Mitchell & Netravali spline")
plt.plot(x, lanczos8(x), color="black",
         linewidth=2.0, linestyle="-", label="Lanczos 8")
plt.plot(x, lanczos10(x), color="gray",
         linewidth=2.0, linestyle="-", label="Lanczos 10")
plt.title("Some kernel functions")
plt.legend()

plt.figure(11)
plt.clf()
ker1 = CardinalCubicSpline(-1/2)
ker2 = CardinalCubicSpline(0)
x = linspace(-3,4,141)
plt.plot(x, ker1(x), color="darkgreen",
         linewidth=2.0, linestyle="-", label="f(x; c = -1/2)")
plt.plot(x, ker1'(x), color="darkblue",
         linewidth=2.0, linestyle="-", label="f'(x; c = -1/2)")
plt.plot(x, ker2(x), color="darkcyan",
         linewidth=2.0, linestyle="-", label="f(x; c = 0)")
plt.plot(x, ker2'(x), color="firebrick",
         linewidth=2.0, linestyle="-", label="f'(x; c = 0)")
plt.title("Cardinal cubic splines")
plt.legend()

plt.figure(12)
plt.clf()
ker = LinearSpline()
x = linspace(-3,4,141)
plt.plot(x, ker(x), color="darkgreen",
         linewidth=2.0, linestyle="-", label="f(x)")
plt.plot(x, ker'(x), color="darkblue",
         linewidth=2.0, linestyle="-", label="f'(x)")
plt.title("Linear B-spline")
plt.legend()

plt.figure(13)
plt.clf()
ker = QuadraticSpline()
x = linspace(-3,4,141)
plt.plot(x, ker(x), color="darkgreen",
         linewidth=2.0, linestyle="-", label="f(x)")
plt.plot(x, ker'(x), color="darkblue",
         linewidth=2.0, linestyle="-", label="f'(x)")
plt.title("Quadratic B-spline")
plt.legend()

plt.figure(14)
plt.clf()
ker = CubicSpline()
x = linspace(-3,4,141)
plt.plot(x, ker(x), color="darkgreen",
         linewidth=2.0, linestyle="-", label="f(x)")
plt.plot(x, ker'(x), color="darkblue",
         linewidth=2.0, linestyle="-", label="f'(x)")
plt.title("Cubic B-spline")
plt.legend()

end # module


text:
"""
An extensible code selection interface.

The `Selectors` module provides an extensible way to write code that has to dispatch on
different predicates without hardcoding the control flow into a single chain of `if`
statements.

In the following example a selector for a simple condition is implemented and the generated
selector code is described:

```julia
abstract type MySelector <: Selectors.AbstractSelector end

# The different cases we want to test.
abstract type One    <: MySelector end
abstract type NotOne <: MySelector end

# The order in which to test the cases.
Selectors.order(::Type{One})    = 0.0
Selectors.order(::Type{NotOne}) = 1.0

# The predicate to test against.
Selectors.matcher(::Type{One}, x)    = x === 1
Selectors.matcher(::Type{NotOne}, x) = x !== 1

# What to do when a test is successful.
Selectors.runner(::Type{One}, x)    = println("found one")
Selectors.runner(::Type{NotOne}, x) = println("not found")

# Test our selector with some numbers.
for i in 0:5
    Selectors.dispatch(MySelector, i)
end
```

`Selectors.dispatch(Selector, i)` will behave equivalent to the following:

```julia
function dispatch(::Type{MySelector}, i::Int)
    if matcher(One, i)
        runner(One, i)
    elseif matcher(NotOne, i)
        runner(NotOne, i)
    end
end
```

and further to

```julia
function dispatch(::Type{MySelector}, i::Int)
    if i === 1
        println("found one")
    elseif i !== 1
        println("not found")
    end
end
```

The module provides the following interface for creating selectors:

- [`order`](@ref)
- [`matcher`](@ref)
- [`runner`](@ref)
- [`strict`](@ref)
- [`disable`](@ref)
- [`dispatch`](@ref)

"""
module Selectors

import InteractiveUtils: subtypes

"""
Root selector type. Each user-defined selector must subtype from this, i.e.

```julia
abstract type MySelector <: Selectors.AbstractSelector end

abstract type First  <: MySelector end
abstract type Second <: MySelector end
```
"""
abstract type AbstractSelector end

"""
Define the precedence of each case in a selector, i.e.

```julia
Selectors.order(::Type{First})  = 1.0
Selectors.order(::Type{Second}) = 2.0
```

Note that the return type must be `Float64`. Defining multiple case types to have the same
order will result in undefined behaviour.
"""
function order end

"""
Define the matching test for each case in a selector, i.e.

```julia
Selectors.matcher(::Type{First}, x)  = x == 1
Selectors.matcher(::Type{Second}, x) = true
```

Note that the return type must be `Bool`.

To match against multiple cases use the [`Selectors.strict`](@ref) function.
"""
function matcher end

"""
Define the code that will run when a particular [`Selectors.matcher`](@ref) test returns
`true`, i.e.

```julia
Selectors.runner(::Type{First}, x)  = println("`x` is equal to `1`.")
Selectors.runner(::Type{Second}, x) = println("`x` is not equal to `1`.")
```
"""
function runner end

"""
Define whether a selector case will "fallthrough" or not when successfully matched against.
By default matching is strict and does not fallthrough to subsequent selector cases.

```julia
# Adding a debugging selector case.
abstract type Debug <: MySelector end

# Insert prior to all other cases.
Selectors.order(::Type{Debug}) = 0.0

# Fallthrough to the next case on success.
Selectors.strict(::Type{Debug}) = false

# We always match, regardless of the value of `x`.
Selectors.matcher(::Type{Debug}, x) = true

# Print some debugging info.
Selectors.runner(::Type{Debug}, x) = @show x
```
"""
strict(::Type{T}) where {T <: AbstractSelector} = true

"""
Disable a particular case in a selector so that it is never used.

```julia
Selectors.disable(::Type{Debug}) = true
```
"""
disable(::Type{T}) where {T <: AbstractSelector} = false

"""
Call `Selectors.runner(T, args...)` where `T` is a subtype of
`MySelector` for which `matcher(T, args...)` is `true`.

```julia
Selectors.dispatch(MySelector, args...)
```
"""
function dispatch(::Type{T}, x...) where T <: AbstractSelector
    for t in (sort(subtypes(T); by = order))
        if !disable(t) && matcher(t, x...)
            runner(t, x...)
            strict(t) && return
        end
    end
    runner(T, x...)
end

end


text:
let
    for use_autodiff in (false, true)
        for (name, prob) in Optim.UnconstrainedProblems.examples
            if prob.isdifferentiable
                f_prob = prob.f
                iterations = if name == "Rosenbrock"
                        5000 # Zig-zagging
                    elseif name == "Powell"
                        80000 # Zig-zagging as the problem is (intentionally) ill-conditioned
                    else
                        1000
                end
                res = Optim.optimize(f_prob, prob.initial_x, GradientDescent(),
                                     OptimizationOptions(autodiff = use_autodiff,
                                                         iterations = iterations))
                @assert norm(res.minimum - prob.solutions, Inf) < 1e-2
            end
        end
    end

    function f_gd_1(x)
      (x[1] - 5.0)^2
    end

    function g_gd_1(x, storage)
      storage[1] = 2.0 * (x[1] - 5.0)
    end

    initial_x = [0.0]

    d = DifferentiableFunction(f_gd_1, g_gd_1)

    results = Optim.optimize(d, initial_x, method=GradientDescent())
    @test_throws ErrorException Optim.x_trace(results)
    @assert Optim.g_converged(results)
    @assert norm(Optim.minimizer(results) - [5.0]) < 0.01

    eta = 0.9

    function f_gd_2(x)
      (1.0 / 2.0) * (x[1]^2 + eta * x[2]^2)
    end

    function g_gd_2(x, storage)
      storage[1] = x[1]
      storage[2] = eta * x[2]
    end

    d = DifferentiableFunction(f_gd_2, g_gd_2)

    results = Optim.optimize(d, [1.0, 1.0], method=GradientDescent())
    @test_throws ErrorException Optim.x_trace(results)
    @assert Optim.g_converged(results)
    @assert norm(Optim.minimizer(results) - [0.0, 0.0]) < 0.01
end


text:
###############################################################################
##### Fit method for a continuous response
###############################################################################

function StatsBase.sample(bs::RegBartState, bm::BartModel)
  posterior = RegBartPosterior(bm)
  @time for s in 1:bm.opts.S
    drawtrees!(bs, bm)
    drawσ!(bs, bm)
    if bm.hypers.sparse
      draws!(bs, bm)
      drawα!(bs, bm)
      # println(bs.shape)
    end
    if s > bm.opts.nburn
      posterior.sdraws[:,s - bm.opts.nburn] = exp.(bs.s)
      posterior.σdraws[s - bm.opts.nburn] = bs.σ
      posterior.treedraws[s - bm.opts.nburn] = [deepcopy(t.tree) for t in bs.ensemble.trees]
    end
    if s % 100 == 0
      println("MCMC iteration $s complete.")
    end
  end
  posterior
end

function StatsBase.fit(BartModel, X::Matrix{Float64}, y::Vector{Float64}, opts = Opts(); hyperags...)
  bm = BartModel(X, y, opts; hyperags...)
  states = RegBartState(bm)
  init_trees = deepcopy(map(state -> Tree[bt.tree for bt in state.ensemble.trees], states))
  post = pmap(bs -> sample(bs, bm), states)
  println("Processing chains...")
  RegBartChain(
    bm,
    init_trees,
    reshape(reduce(hcat, [chain.sdraws for chain in post]), bm.td.p, bm.opts.ndraw, bm.opts.nchains),
    reshape(reduce(vcat, [chain.treedraws for chain in post]), bm.opts.ndraw, 1, bm.opts.nchains),
    reshape(reduce(vcat, [chain.σdraws for chain in post]), bm.opts.ndraw, 1, bm.opts.nchains)
  )
end

function update(post::RegBartChain, ndraw::Int)
  s = post.bm.opts.ndraw
  new_opts = Opts(ndraw = ndraw, nburn = 0, nchains = post.bm.opts.nchains)
  states = RegBartState[]
  for c in 1:post.bm.opts.nchains
    trees = post.treedraws[s,1,c]
    σ = post.σdraws[s,1,c]
    S = [leafprob(post.bm.td.X, tree) for tree in trees]
    fhats = reduce(hcat, [S[t] * getμ(trees[t]) for t in eachindex(trees)])
    yhat = vec(sum(fhats, dims = 2))
    bt = Vector{BartTree}(undef, post.bm.hypers.m)
    for t in eachindex(trees)
      rt = post.bm.td.y - sum(fhats[:,eachindex(trees) .!= t], dims = 2)
      Ω = inv(transpose(S[t]) * S[t] / σ^2 + I / post.bm.hypers.τ)
      rhat = vec(transpose(S[t]) * rt / σ^2)
      bt[t] = BartTree(trees[t], S[t], SuffStats(size(S[t], 2), Ω, rhat))
    end
    bs = RegBartState(BartEnsemble(bt), post.mdraws[:,s,c], σ, ones(post.bm.td.p) ./ post.bm.td.p)
    push!(states, bs)
  end
  bm = BartModel(post.bm.hypers, new_opts, post.bm.td)
  newdraws = pmap(bs -> sample(bs, bm), states)
  RegBartChain(
    BartModel(bm.hypers, Opts(ndraw = ndraw + post.bm.opts.ndraw, nburn = post.bm.opts.nburn), bm.td),
    post.init_trees,
    hcat(post.sdraws,
      reshape(reduce(hcat, [chain.sdraws for chain in newdraws]), bm.td.p, bm.opts.ndraw, bm.opts.nchains)),
    vcat(post.treedraws,
      reshape(reduce(vcat, [chain.treedraws for chain in newdraws]), bm.opts.ndraw, 1, bm.opts.nchains)),
    vcat(post.σdraws,
      reshape(reduce(vcat, [chain.σdraws for chain in newdraws]), bm.opts.ndraw, 1, bm.opts.nchains))
  )
end

###############################################################################
##### Fit method for a binary response ---> ProbitBART
###############################################################################

function StatsBase.sample(bs::ProbitBartState, bm::BartModel)
  posterior = ProbitBartPosterior(bm)
  @time for s in 1:bm.opts.S
    drawtrees!(bs, bm)
    if bm.hypers.sparse
      draws!(bs, bm)
      drawα!(bs, bm)
      # println(bs.shape)
    end
    if s > bm.opts.nburn
      posterior.sdraws[:,s - bm.opts.nburn] = exp.(bs.s)
      posterior.treedraws[s - bm.opts.nburn] = [deepcopy(t.tree) for t in bs.ensemble.trees]
    end
    if s % 100 == 0
      println("MCMC iteration $s complete.")
    end
  end
  posterior
end

function StatsBase.fit(BartModel, X::Matrix{Float64}, y::Vector{Int}, opts = Opts(); hyperags...)
  bm = BartModel(X, y, opts; hyperags...)
  states = ProbitBartState(bm)
  init_trees = deepcopy(map(state -> Tree[bt.tree for bt in state.ensemble.trees], states))
  post = pmap(bs -> sample(bs, bm), states)
  println("Processing chains...")
  ProbitBartChain(
    bm,
    init_trees,
    reshape(reduce(hcat, [chain.sdraws for chain in post]), bm.td.p, bm.opts.ndraw, bm.opts.nchains),
    reshape(reduce(vcat, [chain.treedraws for chain in post]), bm.opts.ndraw, 1, bm.opts.nchains)
  )
end

function update(post::ProbitBartChain, ndraw::Int)
  s = post.bm.opts.ndraw
  new_opts = Opts(ndraw = ndraw, nburn = 0, nchains = post.bm.opts.nchains)
  states = ProbitBartState[]
  z = map(y -> y == 1 ?
    rand(Truncated(Normal(), 0, Inf)) :
    rand(Truncated(Normal(), -Inf, 0)),
    post.bm.td.y
  )
  for c in 1:post.bm.opts.nchains
    trees = post.treedraws[s,1,c]
    S = [leafprob(post.bm.td.X, tree) for tree in trees]
    fhats = reduce(hcat, [S[t] * getμ(trees[t]) for t in eachindex(trees)])
    yhat = vec(sum(fhats, dims = 2))
    bt = BartEnsemble(Vector{BartTree}(undef, post.bm.hypers.m))
    for t in eachindex(trees)
      rt = z - sum(fhats[:,eachindex(trees) .!= t], dims = 2)
      Ω = inv(transpose(S[t]) * S[t] + I / post.bm.hypers.τ)
      rhat = vec(transpose(S[t]) * rt)
      bt.trees[t] = BartTree(trees[t], S[t], SuffStats(size(S[t], 2), Ω, rhat))
    end
    push!(states, ProbitBartState(bt, yhat, z, 1, ones(post.bm.td.p) ./ post.bm.td.p))
  end
  bm = BartModel(post.bm.hypers, new_opts, post.bm.td)
  newdraws = pmap(bs -> sample(bs, bm), states)
  ProbitBartChain(
    BartModel(bm.hypers, Opts(ndraw = ndraw + post.bm.opts.ndraw, nburn = post.bm.opts.nburn), bm.td),
    post.init_trees,
    hcat(post.sdraws,
      reshape(reduce(hcat, [chain.s for chain in newdraws]), bm.td.p, bm.opts.ndraw, bm.opts.nchains)),
    vcat(post.treedraws,
      reshape(reduce(vcat, [chain.treedraws for chain in newdraws]), bm.opts.ndraw, 1, bm.opts.nchains))
  )
end

###############################################################################
##### Fit method for a survival response ---> SurvBART
###############################################################################

function preprocess(X::Matrix{Float64}, time::Vector{Float64}, event::Vector{Int})
  n, p = size(X)
  ut = sort(unique(time))
  yij = Int[]
  tj = Float64[]
  obs = Int[]
  for i in 1:n
    for t in filter(t -> t <= time[i], ut)
      t < time[i] ? push!(yij, 0) : push!(yij, event[i])
      push!(tj, t)
      push!(obs, i)
    end
  end
  ns = counts(obs)
  Xrep = [transpose(reshape(repeat(X[i,:],  ns[i]), p, ns[i])) for i in 1:n]
  Xrep = reduce(vcat, Xrep)
  (yij, hcat(tj, Xrep))
end

function StatsBase.fit(BartModel, X::Matrix{Float64}, time::Vector{Float64}, event::Vector{Int},
  opts = Opts(); hyperags...)
  y, X = preprocess(X, time, event)
  bm = BartModel(X, y, opts; hyperags...)
  states = ProbitBartState(bm)
  init_trees = deepcopy(map(state -> Tree[bt.tree for bt in state.ensemble.trees], states))
  post = pmap(bs -> sample(bs, bm), states)
  println("Processing chains...")
  treedraws = reduce(hcat, pmap(p -> p.treedraws, post))
  pdraws = reduce(hcat, pmap(t -> cdf.(Normal(), predict(t, bm.td.X)), treedraws))
  ut = sort(unique(time))
  obs = [i for i in 1:length(time) for t in filter(t -> t <= time[i], ut)]
  ns = counts(obs)
  indices = [findall(obs .== i) for i in 1:bm.td.n]
  sdraws = [cumprod(1 .- pdraws[idx,:], dims = 1) for idx in indices]
  sdraws = reduce(vcat, sdraws)
  SurvBartChain(
    bm,
    time,
    event,
    init_trees,
    reshape(reduce(hcat, [chain.sdraws for chain in post]), bm.td.p, bm.opts.ndraw, bm.opts.nchains),
    reshape(reduce(vcat, [chain.treedraws for chain in post]), bm.opts.ndraw, 1, bm.opts.nchains),
    sdraws
  )
end


text:
#
# This file is part of the YAActL.jl Julia package, MIT license
#
# Paul Bayer, 2020
#

"""
    _ACT()

Internal actor status variable.

# Fields

1. `dsp::Dispatch`: dispatch mode,
2. `sta::Tuple`: the actor's status variable,
3. `res::Tuple`: the result of the last behavior execution,
4. `bhv::Func` : the behavior function and its internal arguments,
5. `init::Func`: the init function and its arguments,
6. `term::Func`: the terminate function and its arguments,
7. `self::Link`: the actor's (local or remote) self,
8. `name::Symbol`: the actor's registered name.

see also: [`Dispatch`](@ref), [`Func`](@ref), [`Link`](@ref)
"""
mutable struct _ACT
    dsp::Dispatch
    sta::Any
    res::Any
    bhv::Func
    init::Union{Nothing,Func}
    term::Union{Nothing,Func}
    self::Union{Nothing,Link}
    name::Union{Nothing,Symbol}

    _ACT() = new(full, nothing, nothing, Func(),
                 nothing, nothing, nothing, nothing)
end

_terminate!(A::_ACT, code) = !isnothing(A.term) && A.term.f((A.term.args..., code)...; kwargs...)

# actor dispatch on messages
_act(A::_ACT, msg::Set)    = _act(A, msg, msg.x)
_act(A::_ACT, msg::Become) = A.bhv = msg.x
_act(A::_ACT, msg::Exec)   = send!(msg.from, Response(msg.func.f(msg.func.args...; msg.func.kwargs...), A.self))
_act(A::_ACT, msg::Update) = _act(A, msg, Val(msg.s))
_act(A::_ACT, msg::Query)  = _act(A, msg, Val(msg.x))
_act(A::_ACT, msg::Diag)   = send!(msg.from, Response(A, A.self))
function _act(A::_ACT, msg::Init)
    A.init = msg.x
    A.sta  = A.init.f(A.init.args...; A.init.kwargs...)
end
_act(A::_ACT, msg::Stop) = _terminate!(A, msg.x)
_act(A::_ACT, msg::M) where M<:Message = _act(A, Val(A.dsp), msg)

# dispatch on Query message
_act(A::_ACT, msg::Query, ::Val{:sta}) = send!(msg.from, Response(A.sta, A.self))
_act(A::_ACT, msg::Query, ::Val{:res}) = send!(msg.from, Response(A.res, A.self))
_act(A::_ACT, msg::Query, ::Val{:bhv}) = send!(msg.from, Response(A.bhv.f, A.self))
_act(A::_ACT, msg::Query, ::Val{:dsp}) = send!(msg.from, Response(A.dsp, A.self))
_act(A::_ACT, msg::Query, x) = send!(msg.from, Response("$x not available", A.self))

# dispatch on Update message
_act(A::_ACT, msg::Update, ::Val{:sta})  = A.sta  = msg.x
_act(A::_ACT, msg::Update, ::Val{:dsp})  = A.dsp  = msg.x
_act(A::_ACT, msg::Update, ::Val{:self}) = A.self = msg.x
_act(A::_ACT, msg::Update, ::Val{:name}) = A.name = msg.x
_act(A::_ACT, msg::Update, ::Val{:arg})  = A.bhv =
    Func(A.bhv.f, msg.x.args...;
         pairs((; merge(A.bhv.kwargs, msg.x.kwargs)...))...)
_act(A::_ACT, msg::Update, x) = nothing

_tuple(x) = applicable(length, x) ? Tuple(x) : (x,)

# dispatch on Call message
function _act(A::_ACT, ::Val{full}, msg::Call)
    res = A.bhv.f((A.bhv.args..., msg.x...)...; A.bhv.kwargs...)
    A.res = res
    send!(msg.from, Response(res, A.self))
end
function _act(A::_ACT, ::Val{state}, msg::Call)
    res = A.bhv.f((A.sta, msg.x...)...; A.bhv.kwargs...)
    A.res = res
    !isnothing(res) && (A.sta = A.res)
    send!(msg.from, Response(res, A.self))
end
# dispatch on Cast message
function _act(A::_ACT, ::Val{full},  msg::Cast)
    res = A.bhv.f((A.bhv.args..., msg.x...)...; A.bhv.kwargs...)
    A.res = res
end
function _act(A::_ACT, ::Val{state}, msg::Cast)
    res = A.bhv.f((A.sta, msg.x...)...; A.bhv.kwargs...)
    A.res = res
    !isnothing(res) && (A.sta = A.res)
end
# dispatch on other user defined messages
function _act(A::_ACT, ::Val{full}, msg::M) where M<:Message
    res = A.bhv.f((A.bhv.args..., msg)...; A.bhv.kwargs...)
    A.res = res
end
function _act(A::_ACT, ::Val{state}, msg::M) where M<:Message
    res = A.bhv.f((A.sta, msg)...; A.bhv.kwargs...)
    A.res = res
    !isnothing(res) && (A.sta = A.res)
end

# this is the actor loop
function _act(ch::Channel{Message})
    A = _ACT()
    task_local_storage("_ACT", A)
    while true
        msg = take!(ch)
        _act(A, msg)
        msg isa Stop && break
    end
    isnothing(A.name) || call!(_REG, unregister, A.name)
end

"""
```
Actor([lp::LinkParams], bhv::Function, args1...; kwargs...)
Actor(pid::Int, bhv::Function, args1...; kwargs...)
```
Create a new actor. Start a task listening to messages `msg`
sent over the returned self and executing `bhv(args1..., msg; kwargs...)`
for each message. The actor stops if sent [`Stop()`](@ref).

# Arguments
- `[lp::LinkParams]`: parameters for creating the actor,
- `pid::Int`: process `pid` to create the actor on, this can
    also be given with `lp`,
- `bhv`: a function implementing the actor's behavior,
- `args1...`: first arguments to `bhv` (without possible `msg` arguments),
- `kwargs...`: keyword arguments to `bhv`.

# Returns
- a [`Link`](@ref) to a created actor.

see also: [`LinkParams`](@ref)

# Example

```julia
julia> using YAActL, .Threads

julia> act1 = Actor(threadid)               # start an actor who gives its threadid
Channel{Message}(sz_max:32,sz_curr:0)

julia> call!(act1)                          # call it
1
```
"""
function Actor(lp::LinkParams, bhv::F, args::Vararg{Any, N}; kwargs...) where {F<:Function,N}
    lk = lp.pid == myid() ?
        _link(_act, size=lp.size, taskref=lp.taskref, spawn=lp.spawn) :
        _link(lp.pid, _act, size=lp.size, taskref=lp.taskref, spawn=lp.spawn)
    put!(lk.chn, Update(:self, lk))
    become!(lk, bhv, args...; kwargs...)
    return lk
end
Actor(bhv::F, args...; kwargs...) where {F<:Function} =
    Actor(LinkParams(), bhv, args..., kwargs...)
Actor(pid::Int, bhv::F, args...; kwargs...) where {F<:Function} =
    Actor(LinkParams(pid), bhv, args..., kwargs...)

"""
    become(bhv::Function, args...; kwargs...)

Cause your actor to take on a new behavior. This can only be
called from inside an actor/behavior.

# Arguments
- `bhv::Function`: function implementing the new behavior,
- `args...`: arguments to `bhv` (without `msg`),
- `kwargs...`: keyword arguments to `bhv`.
"""
function become(bhv::F, args::Vararg{Any, N}; kwargs...) where {F<:Function,N}
    act = task_local_storage("_ACT")
    act.bhv = Func(bhv, args...; kwargs...)
end

"""
    stop(code=0)

Cause your actor to exit with `code`.
"""
stop(code=0) = send!(_self(), Stop(code))


text:
module cme257hw6

function jset{T <: Number}(R::Function, x::Array{T,1}, y::Array{T,1}, n_iter::Int64 = 100, escape_tol::Float64 = 100.)
    A = zeros(length(x), length(y))
    for i in 1:length(x)
        for j in 1:length(y)
            z = x[i] + im * y[j]
            k = 0
            while k <= n_iter && abs(z) < escape_tol
                z = R(z)
                k = k + 1
            end
            A[i,j] = k
        end
    end
    return A
end

function jset_viz{T<:Number}(R::Function, x::Array{T,1}, y::Array{T,1}, n_iter::Int64, escape_tol::Float64)
    return spy(jset(R::Function, x::Array{T,1}, y::Array{T,1}, n_iter::Int64, escape_tol::Float64))
end

end # module


text:
using titanic
using Test
using titanic: WeakClassifier, Labels, Dataset, classificationerror, out, computeE, gradE, l2norm, linearkernel, polynomialkernel, gaussiankernel

@testset "titanic.jl" begin

    @testset "adaboost" begin
        x1 = [1, 2, 3, 4]
        x2 = [5, 6, 7, 8]
        h1 = WeakClassifier(1, 3, true)
        h2 = WeakClassifier(2, 2, true)
        h3 = WeakClassifier(3, 4, false)
        h4 = WeakClassifier(4, 8, false)

        @test h1(x1) == -1
        @test h1(x2) == 1

        @test h2(x1) == 1
        @test h2(x2) == 1

        @test h3(x1) == 1
        @test h3(x2) == -1

        @test h4(x1) == 1
        @test h4(x2) == -1
    end

    @testset "dataset" begin
        prediction = Labels([1, -1, 1, 1, -1, 1, 1, 1, -1, -1])
        truth = Labels([-1, -1, -1, 1, 1, 1, -1, 1, 1, -1])

        @test classificationerror(truth, prediction) == 0.5

        @test all(out(truth)[2] .== [0, 0, 0, 1, 1, 1, 0, 1, 1, 0])
    end

    @testset "logistic regression" begin
        X = [1.0 1.0 1.0; 1.0 2.0 3.0]
        k = [1, -1, -1]
        w = [1.5, -0.5]

        @test isapprox(computeE(X, k, w), 0.66016; atol=1e-5)

        @test all(isapprox.(gradE(X, k, w), [0.28450, 0.82532]; atol=1e-5))

        @testset "||$(x)|| = √2" for x in [[1, 1], [-1, -1], [-1, 1], [1, -1]]
            @test isapprox(l2norm(x), sqrt(2); atol=1e-5)
        end

        @test all(classify(Dataset(X), [2.0, -1.0])() .== [1, 1, -1])

    end

    @testset "svm" begin
        x1 = [1, 2, 3, 4, 5]
        x2 = [1, 1, 1, 1, 1]

        @test linearkernel(x1, x2) == 15
        @test polynomialkernel(x1, x2; degree=1) == 16
        @test polynomialkernel(x1, x2; degree=2) == 256
        @test isapprox(gaussiankernel(x1, x2; variance=1.0), 3.1e-7; atol=1e-8)
    end

end


text:
## error handling

function tradesyserror(statussignal::Bool)
  if !statussignal
    println("Trading system error encountered")
  end
  return statussignal
end


text:
"""
A form batch is a vectorized form with n primitives transformed into a simpler
representation: one primitive repositioned n times.

On certain backends this leads to more efficient drawing. For example, SVG can
be shortened by using <def> and <use> tags, and raster graphics can render the
form primitive to a back buffer and blit it into place for faster drawing.

Batching is an optimization transform that happens at draw time. There's
currently no mechanism to manually batch. E.g. contexts cannot have FormBatch
children.
"""
struct FormBatch{P <: FormPrimitive}
    primitive::P
    offsets::Vector{AbsoluteVec2}
end

"""
Attempt to batch a form. Return a Nullable{FormBatch} which is null if the Form
could not be batched, and non-null if the original form can be replaced with teh
resulting FormBatch.
"""
batch(form::Form{P}) where P = Nullable{FormBatch{P}}()

# Note: in tests using random data, this optimization wasn't worth it. I'm
# keeping it around out of hopes I find a more clever version that is
# worthwhile, or benchmarks using real data show different results.

# maximum distance between offsets to be considered redundand in mm
const offset_redundancy_threshold = 0.05

"""
Produce a new array of offsets in which near duplicate values have been removed.
"""
function filter_redundant_offsets!(offsets::Vector{AbsoluteVec2})
    isempty(offsets) && return offsets

    sort!(offsets)
    nonredundant_offsets = AbsoluteVec2[offsets[1]]
    for i in 2:length(offsets)
        # use l1 distance for perf
        d = abs(offsets[i-1][1].value - offsets[i][1].value) +
            abs(offsets[i-1][2].value - offsets[i][2].value)
        d > offset_redundancy_threshold && push!(nonredundant_offsets, offsets[i])
    end
    @show (length(offsets), length(nonredundant_offsets))

    return nonredundant_offsets
end

#=
function batch{T <: CirclePrimitive}(form::Form{T})
    # circles can be batched if they all have the same radius.
    r = form.primitives[1].radius
    n = length(form.primitives)
    for i in 2:n
        form.primitives[i].radius == r || return Nullable{FormBatch{CirclePrimitive}}()
    end

    prim = CirclePrimitive((0mm, 0mm), r)
    offsets = Array{AbsoluteVec2}(n)
    for i in 1:n
        offsets[i] = form.primitives[i].center
    end

    return Nullable(FormBatch(prim, offsets))
end
=#

# TODO: same for polygon, rectangle, ellipse

# TODO: batch needs to be exposed as something that users can construct and
# insert into a context. It doesn't make sense to make the same polygon over and
# over and then try to convert it to FormBach.

# Don't attempt to optimize for batching if the form is smaller than this.
const batch_length_threshold = 100

"""
Count the number of unique primitives in a property, stopping when max_count is
exceeded.
"""
function count_unique_primitives(property::Property, max_count::Int)
    unique_primitives = Set{eltype(property.primitives)}()
    for primitive in property.primitives
        push!(unique_primitives, primitive)
        length(unique_primitives) > max_count && break
    end

    return length(unique_primitives)
end

"""
Remove and return vector forms and vector properties from the Context.
"""
function excise_vector_children!(ctx::Context)
    # excise vector forms
    prev_form_child = form_child = ctx.form_children
    forms = Form[]
    while !isa(form_child, ListNull)
        if length(form_child.head.primitives) > 1
            push!(forms, form_child.head)
            if prev_form_child == form_child
                prev_form_child = ctx.form_children = form_child.tail
            else
                prev_form_child.tail = form_child.tail
            end
        else
            prev_form_child = form_child
        end

        form_child = form_child.tail
    end

    # excise vector properties
    prev_property_child = property_child = ctx.property_children
    properties = Property[]
    while !isa(property_child, ListNull)
        if length(property_child.head.primitives) > 1
            push!(properties, property_child.head)
            if prev_property_child == property_child
                prev_property_child = ctx.property_children = property_child.tail
            else
                prev_property_child.tail = property_child.tail
            end
        else
            prev_property_child = property_child
        end

        property_child = property_child.tail
    end

    return (forms, properties)
end

"""
Attempt to transform a tree into an equivalent tree that can more easily be
batched.

What this does is look for patterns in which a long vector form is accompanied
by a large vector property that has a relatively small number of unique values.
If there are n unique values, we can split it into n contexts, each with a
shorter vector form and only scalar properties.
"""
function optimize_batching(ctx::Context)
    # condition 1: has a 1 or more long vector forms
    max_form_length = 0
    form_child = ctx.form_children
    while !isa(form_child, ListNull)
        max_form_length = max(max_form_length, length(form_child.head.primitives))
        form_child = form_child.tail
    end

    max_form_length < batch_length_threshold && return ctx

    # condition 2: has a 1 or more long vector properties each with a smaller
    # number of unique values
    max_count = div(max_form_length, batch_length_threshold) + 1
    max_unique_primitives = 0
    prop_child = ctx.property_children
    while !isa(prop_child, ListNull)
        if length(prop_child.head.primitives) > 1
            max_unique_primitives = max(max_unique_primitives,
                    count_unique_primitives(prop_child.head, max_count))
        end
        prop_child = prop_child.tail
    end

    # don't batch when there are not many forms per unique property primitive
    if max_unique_primitives == 0 ||
        div(max_form_length, max_unique_primitives) + 1 < batch_length_threshold
        return ctx
    end

    # non-destructive since this happens at draw time and draw should not modify
    # the context.
    ctx = copy(ctx)

    # step 1: remove vector form and vector properties
    forms, properties = excise_vector_children!(ctx)

    # step 2: split primitives into groups on the cross product of property
    # primives
    n = length(forms[1].primitives)
    grouped_forms = Dict{UInt, Vector{Form}}()
    grouped_properties = Dict{UInt, Vector{Property}}()
    for i in 1:n
        h = UInt(0)
        for property in properties
            h = hash(property.primitives[i], h)
        end

        if !haskey(grouped_forms, h)
            grouped_forms[h] = Form[similar(form) for form in forms]
            group_prop = Array{Property}(length(properties))
            for j in 1:length(properties)
                group_prop[j] = Property([properties[j].primitives[i]])
            end
            grouped_properties[h] = group_prop
        end

        for j in 1:length(forms)
            push!(grouped_forms[h][j].primitives, forms[j].primitives[i])
        end
    end

    # step 3: put forms in new contexts and insert into the ctx
    for (h, fs) in grouped_forms
        subctx = context()
        compose!(subctx, fs...)
        compose!(subctx, grouped_properties[h]...)
        compose!(ctx, subctx)
    end

    return ctx
end


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

name = "libint"
version = v"2.7.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/evaleev/libint/releases/download/v2.7.0-beta.5/libint-2.7.0-beta.5.tgz", "38c630d9b4433340079ebc7a93f05d7a285833883327bfc35dd208249e06f05f"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libint-*/
mkdir build && cd build/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DLIBINT2_BUILD_SHARED_AND_STATIC_LIBS=ON ..
cmake --build . --target install -- -j${nproc}
# Remove large static archive
rm ${prefix}/lib/liblibint2-static.a
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter!(p -> !isa(p, Windows), supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("liblibint2", :Libint2)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency("Eigen_jll"),
    BuildDependency("boost_jll"),
]

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


text:
"""
    fit(::Type{<:AbstractModel}, df; kwargs)

Optimize the parameters of a model using measurements in `df` and the initialisation values in
`kwargs`. Note that the columns in `df` should match exactly the names and units used in the
model. See particular implementations for more details.
"""
function fit end

# See other files for the implementation of specific model fitting


text:
struct Travelers end

n_agents(simpleGame::Travelers) = 2

ordered_actions(simpleGame::Travelers, i::Int) = 2:100
ordered_joint_actions(simpleGame::Travelers) = vec(collect(Iterators.product([ordered_actions(simpleGame, i) for i in 1:n_agents(simpleGame)]...)))

n_joint_actions(simpleGame::Travelers) = length(ordered_joint_actions(simpleGame))
n_actions(simpleGame::Travelers, i::Int) = length(ordered_actions(simpleGame, i))

function reward(simpleGame::Travelers, i::Int, a)
    if i == 1
        noti = 2
    else
        noti = 1
    end
    if a[i] == a[noti]
        r = a[i]
    elseif a[i] < a[noti]
        r = a[i] + 2
    else
        r = a[noti] - 1
    end
    return r
end

function joint_reward(simpleGame::Travelers, a)
    return [reward(simpleGame, i, a) for i in 1:n_agents(simpleGame)]
end

function SimpleGame(simpleGame::Travelers)
    return SimpleGame(
        0.9,
        vec(collect(1:n_agents(simpleGame))),
        [ordered_actions(simpleGame, i) for i in 1:n_agents(simpleGame)],
        (a) -> joint_reward(simpleGame, a)
    )
end


text:
module GPUArrays

using Serialization
using Random
using LinearAlgebra
using Printf

using LinearAlgebra.BLAS
using Base.Cartesian

using AbstractFFTs

using Adapt

# device array
include("device/abstractarray.jl")
include("device/indexing.jl")
include("device/synchronization.jl")

# host array
include("host/abstractarray.jl")
include("host/devices.jl")
include("host/execution.jl")
include("host/construction.jl")
## integrations and specialized functionality
include("host/base.jl")
include("host/indexing.jl")
include("host/broadcast.jl")
include("host/mapreduce.jl")
include("host/linalg.jl")
include("host/random.jl")
include("host/quirks.jl")

# CPU reference implementation
include("array.jl")

end # module


text:
using EntanglementDist

# Change the following to perform further tests:
eps = 0.2;

# rho = wernerState(0.9);
# rho = rState(0.7);

# Customized state:
p = 0.7
rho = sState(p);

# Output initial fidelity
Finit = entFidelity(rho);
print("Initial fidelity is F=", round(Finit,3), "\n");

# Test the seesaw optimization for the filtering scheme
# Adapted to work well for states of the form p EPR + (1-p) |01><01|
print("Testing the seesaw method\n");
print("Epsilon is set to ", eps, "\n");

print("Computing Choi states\n");
(CA,CB) = filteringMakeChoi(eps);
print("Compute the filtering performance without Choi state to double check:\n");
(rhoQC, P, Fwo, pwo) = filtering(rho, eps);
print("Checked F=",round(Fwo,3), " and psucc=",round(pwo,3),"\n");

print("Run the SEESAW Optimization\n");
(newCA, newCB, Fnew, pnew) = seesaw(rho,2,2,2,CA, CB, P);