text:
module SpiderMonkey

export smo

using Statistics, Distributions, Random

"""
smo(D::Int,P::Int,f::Function,ɳ = 0.5,λ = 25,α = 20,β = 8,τ = 100,ub = fill(100,(D,)),lb = fill(-100,(D,)),verbose = false)

Returns a tuple of x_opt (Array of values at which function is optimum) and f_opt (optimum function value)

Spider Monkey Optimization (SMO) is a global optimization algorithm inspired by Fission-Fusion social (FFS) structure of spider monkeys during their foraging behavior. SMO exquisitely depicts two fundamental concepts of swarm intelligence: self-organization and division of labor.

# Arguments
- `D`: Dimension of the function being optimized, must be an integer.
- `P`: Population, Number of Spider Monkeys for the heuristic search, must be an integer.
- `f`: function to be optimized, example reference - BlackBoxOptimizationBenchmarking functions.
- `ɳ`: perturbation rate
- `λ`: local leader limit, max count for which the position of the local leader remains unchanged
- `α`: global leader limit, max count for which the position of the global leader remains unchanged
- `β`: maximum number of groups allowed
- `τ`: maximum number of iterations
- `ub`: upper bound values for spider monkey position coordinates, vector of length D.
- `lb`: lower bound values for spider monkey position coordinates, vector of length D.
- `verbose`: print to console the iteration number, Global leader fitness (f_opt) and number of groups for each iteration.
# Examples
## Optimize a function
Find out the optimum value of a function (numeric optimization) and the value as a vector of dimention D at which the function has optimum value.
```julia-repl
julia> using BlackBoxOptimizationBenchmarking;
julia> ub = fill(5,(100,));
julia> lb = fill(-5,(100,));
julia> smo(100,30,BlackBoxOptimizationBenchmarking.F1.f,0.5,25,20,8,100,ub,lb);
```
"""

mutable struct SMGroup{First,Population,LocalLeader,LLCount,LLIndex}
    first::First
    last::Population
    local_leader::LocalLeader
    local_leader_count::LLCount
    local_leader_index::LLIndex
    SMGroup(First,Population,LocalLeader,LLCount,LLIndex) = new{typeof(First),typeof(Population),typeof(LocalLeader),typeof(LLCount),typeof(LLIndex)}(First,Population,LocalLeader,LLCount,LLIndex)
end

function smo(
        D::Int, # Dimensions
        P::Int, # Population
        f::Function, # Function to optimize (in example format)
        ɳ = 0.5, # Perturbation rate
        λ = 25, # Local leader Limit
        α = 20, # Global Leader Limit
        β = 8, # Max Group Count
        τ = 100, # Iterations
        ub = fill(100,(D,)), # Upper bound for value
        lb = fill(-100,(D,)), # Lower bound for values
        verbose = false
        )

        groups = []
        GLC = 0
        # initialization of monkey population
        SM = rand(P,D)
        for i in 1:P
            for j in 1:D
                SM[i,j] = lb[j] + rand(Uniform(),1)[1]*(ub[j] - lb[j])
            end
        end

        fitness = map(f,eachrow(SM))
        LLIndex = argmin(fitness)
        # Local and Global leaders same with one group
        LL = GL = SM[LLIndex,:]
        group = SMGroup(1,P,LL,0,LLIndex)
        # Global Leader fitness
        GL_fitness = minimum(fitness)
        push!(groups,group)

        while τ > 0
            if verbose
                println("Iteration: $(τ), Global Leader Fitness: $(GL_fitness), Groups: $(length(groups))")
            end
            SM_temp = copy(SM)
            prev_fitness = map(f,eachrow(SM_temp))
            SM = local_leader_phase(groups,D,ɳ,SM,lb,ub)
            new_fitness = map(f,eachrow(SM))
            for i in 1:P
                if new_fitness[i] > prev_fitness[i]
                    SM[i,:] = SM_temp[i,:]
                    fitness[i] = prev_fitness[i]
                else
                    fitness[i] = new_fitness[i]
                end
            end
            prev_fitness = copy(fitness)
            SM_temp = copy(SM)
            SM = global_leader_phase(groups,D,GL,fitness,SM,lb,ub)
            new_fitness = map(f,eachrow(SM))
            for i in 1:P
                if new_fitness[i] > prev_fitness[i]
                    SM[i,:] = copy(SM_temp[i,:])
                    fitness[i] = prev_fitness[i]
                else
                    fitness[i] = new_fitness[i]
                end
            end
            for group in groups
                temp_LL = copy(group.local_leader)
                group.local_leader = SM[argmin(fitness),:]
                if all(temp_LL .== group.local_leader)
                    group.local_leader_count += 1
                end
            end
            temp_GL = copy(GL)
            glb_check_SM = SM
            GL = SM[argmin(fitness),:]
            GL_fitness = minimum(fitness)
            if all(temp_GL .== GL)
                GLC += 1
            end
            SM, groups = local_leader_decision(groups,D,ɳ,λ,SM,lb,ub,GL)
            groups = global_leader_decision(groups,D,GLC,fitness,SM,lb,ub,β,α,P)
            τ -= 1
        end
        return GL,GL_fitness
    end

function local_leader_decision(groups,dimensions,ɳ,λ,SM,lb,ub,GL)
    for group in groups
        if group.local_leader_count > λ
            group.local_leader_count = 0

            GS = group.last - group.first + 1
	        for i in 1:GS
	            for j in 1:dimensions
	                if rand(Uniform(),1)[1] >= ɳ
	                    SM[i,j] += lb[j] + rand(Uniform(),1)[1]*(ub[j] - lb[j])
	                else
	                    SM[i,j] += rand(Uniform(),1)[1]*(GL[j]-SM[i,j])+rand(Uniform(),1)[1]*(SM[i,j]-group.local_leader[j])
                    end
	                if SM[i,j] > ub[j]
	                    SM[i,j] = ub[j]
	                elseif SM[i,j] < lb[j]
	                    SM[i,j] = lb[j]
                    end
                end
            end
        end
    end
    return SM,groups
end

function global_leader_decision(groups,dimensions,GLC,fitness,SM,lb,ub,β,α,P)
    if GLC > α
        GLC = 0
        if length(groups) < β
            g1 = rand(groups)
            if g1.last > g1.first
                temp_last = g1.last
                g1.last = g1.first + cld((g1.last - g1.first),2) - 1
                g1.local_leader_count = 0
                g2 = SMGroup(g1.last + 1,temp_last,SM[argmin(fitness[g1.last:temp_last]),:],0,argmin(fitness[g1.last:temp_last]))
                g1.local_leader = SM[argmin(fitness[g1.first:g1.last]),:]
                push!(groups,g2)
            end
        else
            deleteat!(groups,[i for i in 1:length(groups)])
            g = SMGroup(1,P,SM[argmin(fitness),:],0,argmin(fitness))
            push!(groups,g)
        end
    end
    return groups
end

function global_leader_phase(groups,dimensions,GL,fitness,SM,lb,ub)
    for group in groups
        count = 1
        GS = group.last - group.first + 1
        while count < GS
            for i in 1:GS
                if rand(Uniform(),1)[1] < prob(fitness,group,i)
                    count += 1
                    j = rand(1:dimensions)
                    r = rand(group.first:group.last)
                    if r == i
                        r = i == 1 ? r+1 : r-1
                    end
                    SM[i,j] += rand(Uniform(),1)[1]*(GL[j]-SM[i,j]) + rand(Uniform(-1,1),1)[1]*(SM[r,j]-SM[i,j])
                    if SM[i,j] > ub[j]
                        SM[i,j] = ub[j]
                    elseif SM[i,j] < lb[j]
                    	SM[i,j] = lb[j]
                    end
                end
            end
        end
    end
    return SM
end

function local_leader_phase(groups,dimensions,ɳ,SM,lb,ub)

    for group in groups
        for i in group.first:group.last
            for j in 1:dimensions
                if rand(Uniform(),1)[1] >= ɳ
                    r = rand(group.first:group.last)
                    if r == i
                        r = i == 1 ? r+1 : r-1
                    end
                    SM[i,j] += rand(Uniform(),1)[1]*(group.local_leader[j] - SM[i,j]) + rand(Uniform(-1,1),1)[1]*(SM[r,j] - SM[i,j])
                    if SM[i,j] > ub[j]
                        SM[i,j] = ub[j]
                    elseif SM[i,j] < lb[j]
                        SM[i,j] = lb[j]
                    end
                end
            end
        end
    end
    return SM
end

function prob(fitness,group,i)
    return 0.9*(fitness[group.local_leader_index]/fitness[i]) + 0.1
end

end # module


text:
export BasicBlock, unsafe_delete!,
       terminator, name,
       move_before, move_after

@checked struct BasicBlock <: Value
    ref::API.LLVMValueRef
end
register(BasicBlock, API.LLVMBasicBlockValueKind)

BasicBlock(ref::API.LLVMBasicBlockRef) = BasicBlock(API.LLVMBasicBlockAsValue(ref))
Base.unsafe_convert(::Type{API.LLVMBasicBlockRef}, bb::BasicBlock) = API.LLVMValueAsBasicBlock(bb)

# forward declarations
@checked struct Function <: GlobalObject
    ref::API.LLVMValueRef
end

BasicBlock(f::Function, name::String; ctx::Context) =
    BasicBlock(API.LLVMAppendBasicBlockInContext(ctx, f, name))
BasicBlock(bb::BasicBlock, name::String; ctx::Context) =
    BasicBlock(API.LLVMInsertBasicBlockInContext(ctx, bb, name))

unsafe_delete!(::Function, bb::BasicBlock) = API.LLVMDeleteBasicBlock(bb)
Base.delete!(::Function, bb::BasicBlock) =
    API.LLVMRemoveBasicBlockFromParent(bb)

parent(bb::BasicBlock) = Function(API.LLVMGetBasicBlockParent(bb))

terminator(bb::BasicBlock) = Instruction(API.LLVMGetBasicBlockTerminator(bb))

name(bb::BasicBlock) =
    unsafe_string(API.LLVMGetBasicBlockName(bb))

move_before(bb::BasicBlock, pos::BasicBlock) =
    API.LLVMMoveBasicBlockBefore(bb, pos)
move_after(bb::BasicBlock, pos::BasicBlock) =
    API.LLVMMoveBasicBlockAfter(bb, pos)

## instruction iteration

export instructions

struct BasicBlockInstructionSet
    bb::BasicBlock
end

instructions(bb::BasicBlock) = BasicBlockInstructionSet(bb)

Base.eltype(::BasicBlockInstructionSet) = Instruction

function Base.iterate(iter::BasicBlockInstructionSet,
                      state=API.LLVMGetFirstInstruction(iter.bb))
    state == C_NULL ? nothing : (Instruction(state), API.LLVMGetNextInstruction(state))
end

Base.last(iter::BasicBlockInstructionSet) =
    Instruction(API.LLVMGetLastInstruction(iter.bb))

Base.isempty(iter::BasicBlockInstructionSet) =
    API.LLVMGetLastInstruction(iter.bb) == C_NULL

Base.IteratorSize(::BasicBlockInstructionSet) = Base.SizeUnknown()

ext: jl


text:
# Second-Order upwind central scheme
# Based on:
# Kurganov A., Noelle S., Petrova G., Semidiscrete Central-Upwind schemes
# for hyperbolic Conservation Laws and Hamilton-Jacobi Equations. SIAM. Sci Comput,
# Vol 23, No 3m pp 707-740. 2001

struct FVCUScheme{ltype <: AbstractSlopeLimiter} <: AbstractFVAlgorithm
  slopeLimiter :: ltype
end

function FVCUScheme(;slopeLimiter=GeneralizedMinmodLimiter())
  FVCUScheme(slopeLimiter)
end

function update_flux_value(uold,∇u,j,mesh,Flux, alg::FVCUScheme,nonscalar::Bool)
    uminus=cellval_at_left(j,uold,mesh)+0.5*cellval_at_left(j,∇u,mesh)
    uplus=cellval_at_right(j,uold,mesh)-0.5*cellval_at_right(j,∇u,mesh)

    if nonscalar
      λm = eigvals(Flux.Dflux(uminus))
      λp = eigvals(Flux.Dflux(uplus))
      aa_plus=maximum((maximum(λm), maximum(λp),0))
      aa_minus=minimum((minimum(λm), minimum(λp),0))
    else
      λm = Flux.Dflux(uminus)
      λp = Flux.Dflux(uplus)
      aa_plus=maximum((λm, λp,0))
      aa_minus=minimum((λm, λp,0))
    end
    # Update numerical fluxes
    if abs(aa_plus-aa_minus) < 1e-8
      return 0.5*(Flux(uminus)+Flux(uplus))
    else
      return (aa_plus*Flux(uminus)-aa_minus*Flux(uplus))/(aa_plus-aa_minus) +
      (aa_plus*aa_minus)/(aa_plus-aa_minus)*(uplus - uminus)
    end
end

# Dissipation Reduced Central upwind Scheme: Second-Order
# Based on:
# Kurganov A., Lin C., On the reduction of Numerical Dissipation in Central-Upwind
# Schemes, Commun. Comput. Phys. Vol 2. No. 1, pp 141-163, Feb 2007.

struct FVDRCUScheme{ltype <: AbstractSlopeLimiter} <: AbstractFVAlgorithm
  slopeLimiter :: ltype
end

function FVDRCUScheme(;slopeLimiter=GeneralizedMinmodLimiter())
  FVDRCUScheme(slopeLimiter)
end

function update_flux_value(uold,∇u,j,mesh,Flux, alg::FVDRCUScheme,nonscalar::Bool)
  uminus=cellval_at_left(j,uold,mesh)+0.5*cellval_at_left(j,∇u,mesh)
  uplus=cellval_at_right(j,uold,mesh)-0.5*cellval_at_right(j,∇u,mesh)

  if nonscalar
    λm = eigvals(Flux.Dflux(uminus))
    λp = eigvals(Flux.Dflux(uplus))
    aa_plus=maximum((λm[end], λp[end],0))
    aa_minus=minimum((λm[1], λp[1],0))
  else
    λm = Flux.Dflux(uminus)
    λp = Flux.Dflux(uplus)
    aa_plus=maximum((λm, λp,0))
    aa_minus=minimum((λm, λp,0))
  end

  # Update numerical fluxes
  if abs(aa_plus-aa_minus) < 1e-8
    return 0.5*(Flux(uminus)+Flux(uplus))
  else
    flm = Flux(uminus); flp = Flux(uplus)
    wint = 1/(aa_plus-aa_minus)*(aa_plus*uplus-aa_minus*uminus-(flp-flm))
    qj = minmod.((uplus-wint)/(aa_plus-aa_minus),(wint-uminus)/(aa_plus-aa_minus))
    return (aa_plus*flm-aa_minus*flp)/(aa_plus-aa_minus) +
            (aa_plus*aa_minus)*((uplus - uminus)/(aa_plus-aa_minus) - qj)
  end
end

##################3
FVCUSchemes = Union{FVCUScheme,FVDRCUScheme}
"""
compute_fluxes!(hh, Flux, uold, mesh, dt, M, alg::FVCUScheme, ::Type{Val{true}})
Numerical flux of Second-Order upwind central scheme in 1D
"""
function compute_fluxes!(fluxes, Flux, uold, mesh, dt, alg::T, nonscalar::Bool, ::Type{Val{true}}) where {T <: FVCUSchemes}
  slopeLimiter = alg.slopeLimiter
  #update vector
  # 1. slopes
  ∇u = compute_slopes(uold, mesh, slopeLimiter, nonscalar, Val{true})

  Threads.@threads for j in node_indices(mesh)
      if nonscalar
          fluxes[:,j] .= update_flux_value(uold,∇u,j,mesh,Flux,alg,nonscalar)
      else
          fluxes[j] = update_flux_value(uold,∇u,j,mesh,Flux,alg,nonscalar)
      end
  end
end

function compute_fluxes!(fluxes, Flux, uold, mesh, dt, alg::T, nonscalar::Bool, ::Type{Val{false}}) where {T <: FVCUSchemes}
  slopeLimiter = alg.slopeLimiter
  #update vector
  # 1. slopes
  ∇u = compute_slopes(uold, mesh, slopeLimiter, nonscalar, Val{false})

  for j in node_indices(mesh)
      if nonscalar
          fluxes[:,j] .= update_flux_value(uold,∇u,j,mesh,Flux,alg,nonscalar)
      else
          fluxes[j] = update_flux_value(uold,∇u,j,mesh,Flux,alg,nonscalar)
      end
  end
end

# Dissipation Reduced Central upwind Scheme: Fifth-Order
# Based on:
# Kurganov A., Lin C., On the reduction of Numerical Dissipation in Central-Upwind
# Schemes, Commun. Comput. Phys. Vol 2. No. 1, pp 141-163, Feb 2007.
# Kurganov, Liu, New adaptive artificial viscosity method for hyperbolic systems
# of conservation laws

struct FVDRCU5Scheme{ltype <: AbstractSlopeLimiter} <: AbstractFVAlgorithm
  slopeLimiter :: ltype
end

function FVDRCU5Scheme(;slopeLimiter=GeneralizedMinmodLimiter())
  FVDRCU5Scheme(slopeLimiter)
end

function compute_fluxes!(fluxes, Flux, uold, mesh, dt, alg::FVDRCU5Scheme, nonscalar::Bool, ::Type{Val{true}})
  #update vector
  Threads.@threads for j in node_indices(mesh)
      if nonscalar
          fluxes[:,j] .= update_flux_value(uold,j,mesh,Flux,alg,nonscalar)
      else
          fluxes[j] = update_flux_value(uold,j,mesh,Flux,alg,nonscalar)
      end
  end
end

function compute_fluxes!(fluxes, Flux, uold, mesh, dt, alg::FVDRCU5Scheme, nonscalar::Bool, ::Type{Val{false}})
  #update vector
  for j in node_indices(mesh)
      if nonscalar
          fluxes[:,j] .= update_flux_value(uold,j,mesh,Flux,alg,nonscalar)
      else
          fluxes[j] = update_flux_value(uold,j,mesh,Flux,alg,nonscalar)
      end
  end
end

function update_flux_value(uold,j,mesh,Flux,alg::FVDRCU5Scheme, nonscalar)
    slopeLimiter = alg.slopeLimiter
    # A fifth-order piecewise polynomial reconstruction
    @inbounds ulll=cellval_at_left(j-2,uold,mesh)
    @inbounds ull=cellval_at_left(j-1,uold,mesh)
    @inbounds ul=cellval_at_left(j,uold,mesh)
    @inbounds ur=cellval_at_right(j,uold,mesh)
    @inbounds urr=cellval_at_right(j+1,uold,mesh)
    @inbounds urrr=cellval_at_right(j+2,uold,mesh)

    uminus = 1/60*(2*ulll-13*ull+47*ul+27*ur-3*urr)
    uplus = 1/60*(-3*ull+27*ul+47*ur-13*urr+2*urrr)

    #Remark: some extrange bug wth eigvals force me to use Lapack
    if nonscalar
      λm = sort(LAPACK.geev!('N','N',Array(Flux.Dflux(uminus)))[1])
      λp = sort(LAPACK.geev!('N','N',Array(Flux.Dflux(uplus)))[1])
      aa_plus=maximum((λm[end], λp[end],0))
      aa_minus=minimum((λm[1], λp[1],0))
    else
      λm = Flux.Dflux(uminus)
      λp = Flux.Dflux(uplus)
      aa_plus=maximum((λm, λp,0))
      aa_minus=minimum((λm, λp,0))
    end
    # Update numerical fluxes
    if abs(aa_plus-aa_minus) < 1e-8
      return 0.5*(Flux(uminus)+Flux(uplus))
    else
      flm = Flux(uminus); flp = Flux(uplus)
      wint = 1/(aa_plus-aa_minus)*(aa_plus*uplus-aa_minus*uminus-
      (flp-flm))
      qj = minmod.((uplus-wint)/(aa_plus-aa_minus),(wint-uminus)/(aa_plus-aa_minus))
      return (aa_plus*flm-aa_minus*flp)/(aa_plus-aa_minus) +
      (aa_plus*aa_minus)*((uplus - uminus)/(aa_plus-aa_minus) - qj)
    end
end


text:
module BeeEncoder

import Base: -, +, *, /, mod, max, min, ==, <, <=, !=, >, >=, sum, show, convert, isequal, in

using Suppressor

export BeeInt,
    BeeBool, BeeModel, BeeSolution, GBL_MODEL,

    @beeint, @beebool, beeint, beebool,

    boolnum, intnum, hasvar, hasbool, hasint, fetchbool, getbool, getint,

    render, reset, solve,

    constrain, @constrain, and, or, iff, alldiff,

    capture_render

# -------------------------------------------------------------
# abstract types
# -------------------------------------------------------------

"An object that can be rendered to BEE syntax"
abstract type BeeObject end
"A symbol in BEE sytanx. Can be either a variable or a value."
abstract type BeeSymbol <: BeeObject end
abstract type BeeBoolean <: BeeSymbol end
abstract type BeeInteger <: BeeSymbol end

BB = Union{BeeBoolean, Bool}
ZZ = Union{BeeInteger, Int}

"Reload `isequal` to allow using `BeeSymbol` comparing them for `Dict`."
isequal(v1::BeeSymbol, v2::BeeSymbol) = v1 === v2

"Reload `in` to allow checking if a `BeeSymbol` is in an array"
function in(v::T, vs::Array{T,1}) where T <: BeeSymbol
    for s in vs
        if v === s
            return true
        end
    end
    false
end

# -------------------------------------------------------------
# BEE integer variable
# -------------------------------------------------------------

"An integer variable in BEE syntax"
struct BeeInt <: BeeInteger
    name::String
    lo::Int
    hi::Int
    function BeeInt(model, name, lo, hi)
        if lo > hi
            error("$lo > $hi")
        end
        if hasvar(model, name)
            error("Variable $name has already been defined in $model")
        end
        var = new(name, lo, hi)
        model.intdict[name] = var
    end
end

BeeInt(name::String, lo::Int, hi::Int) = BeeInt(GBL_MODEL, name, lo, hi)
BeeInt(name::Symbol, lo, hi) = BeeInt(string(name), lo, hi)

macro beeint(name, lo, hi)
    q = Expr(:block)
    if isa(name, Symbol)
        ex = :($(esc(name)) = beeint($(string(name)), $lo, $hi))
        push!(q.args, ex)
        push!(q.args, escvar(name)) # return all of the variables we created
    elseif isa(name, Expr) && name.head == :ref
        vhead = name.args[1]
        vlist = []
        for i in eval(name.args[2])
            vhead = name.args[1]
            vstr = "$vhead$i"
            vsym = Symbol(vstr)
            ex = :($(esc(vsym)) = beeint($(vstr), $lo, $hi))
            push!(vlist, vsym)
            push!(q.args, ex)
        end
        push!(q.args, escvar(vlist))
    end
    q
end

"Create an integer varaible called `name` in `GBL_MODEL`"
beeint(name, lo, hi) = beeint(GBL_MODEL, name, lo, hi)

"Create an integer varaible called `name` in `GBL_MODEL`"
beeint(model, name, lo, hi) = BeeInt(model, name, lo, hi)

render(io::IO, var::BeeInt) = print(io, "new_int($var, $(var.lo), $(var.hi))\n")

# -------------------------------------------------------------
# BEE boolean variable
# -------------------------------------------------------------

"A boolean variable in BEE syntax"
struct BeeBool <: BeeBoolean
    name::String
    function BeeBool(model, name)
        if hasvar(model, name)
            error("Variable $name has already been defined in $model")
        end
        var = new(name)
        model.booldict[name] = var
    end
end
BeeBool(name::String) = BeeBool(GBL_MODEL, name)
BeeBool(name::Symbol) = BeeBool(string(name))

function escvar(var)
    if isa(var, Symbol)
        vs = esc(var)
    else
        vs = Expr(:vect, [esc(v) for v in var]...)
    end
    vs
end

macro beebool(namelist...)
    q = Expr(:block)
    varlist = [] # list of boolean variables to create
    for name in namelist
        if isa(name, Symbol)
            push!(varlist, name)
            ex = :($(esc(name)) = beebool($(string(name))))
            push!(q.args, ex)
        elseif isa(name, Expr) && name.head == :ref
            vhead = name.args[1]
            vlist = []
            for i in eval(name.args[2])
                vhead = name.args[1]
                vstr = "$vhead$i"
                vsym = Symbol(vstr)
                ex = :($(esc(vsym)) = beebool($(vstr)))
                push!(vlist, vsym)
                push!(q.args, ex)
            end
            push!(varlist, vlist)
        end
    end
    if length(varlist) > 1
        ret = Expr(:tuple, map(escvar, varlist)...)
    else
        ret = escvar(varlist[1])
    end
    push!(q.args, ret) # return all of the variables we created
    q
end

"Create a boolean varaible called `name` in `GBL_MODEL`"
beebool(name) = beebool(GBL_MODEL, name)

"Create a boolean varaible called `name` in `model`"
beebool(model, name) = BeeBool(model, name)

render(io::IO, var::BeeBool) = print(io, "new_bool($var)\n")

"The negate of a boolean variable in BEE syntax"
struct BeeNegateBool <:BeeBoolean
    boolvar::BeeBool
end

show(io::IO, v::BeeNegateBool) = print(io, "-", v.boolvar)

(-)(var::BeeBool) = BeeNegateBool(var)
(-)(var::BeeNegateBool) = var.boolvar

# -------------------------------------------------------------
# BEE literals
# -------------------------------------------------------------

"An integer value like `1`"
struct BeeIntLiteral <: BeeInteger
    val::Int
end
show(io::IO, v::BeeIntLiteral) = print(io, v.val)
convert(::Type{BeeInteger}, v::Int) = BeeIntLiteral(v)
convert(::Type{BeeSymbol}, v::Int) = BeeIntLiteral(v)

"A boolean value, i.e., `true` or `false`"
struct BeeBoolLiteral <: BeeBoolean
    val::Bool
end
show(io::IO, v::BeeBoolLiteral) = print(io, v.val)
convert(::Type{BeeBoolean}, v::Bool) = BeeBoolLiteral(v)
convert(::Type{BeeSymbol}, v::Bool) = BeeBoolLiteral(v)

"""
    render(obj::BeeObject)

Render `obj` to BEE syntax and print it to `stdout`.
"""
render(obj::BeeObject) = render(Base.stdout, obj)

render(arr::Array{T, 1}) where T <: Any = render(Base.stdout, arr)

render(io::IO, arr::Array{T, 1}) where T <: BeeObject = show(io, arr)

show(io::IO, v::BeeSymbol) = print(io, v.name)

function show(io::IO, arr::Array{T, 1}) where T <: BeeObject
    print(io, "[", join(arr, ", "), "]")
end

# -------------------------------------------------------------
# BEE expressions
# -------------------------------------------------------------

"""
An expression can be made part of a `BeeObject`, but they themselves cannot be rendered.
"""
abstract type BeeExpression end

function show(io::IO,
                   tuple::NTuple{N, T} where {N, T <: BeeSymbol})
    print(io, join(tuple, ", "))
end

# -------------------------------------------------------------
# BEE Constrains
# -------------------------------------------------------------

struct BeeConstraint <: BeeObject
    name::String
    # Don't check the type of the array
    varlist::Array{T, 1} where T
end

"""
    render(io, cons)

Render `cons` to BEE syntax. For constraints, there's is no difference between how they are rendered
and printed.
"""
render(io::IO, constraint::BeeConstraint) = print(io, constraint, "\n")

function show(io::IO, constraint::BeeConstraint)
    varstr = join(constraint.varlist, ", ")
    print(io, constraint.name, "(", varstr, ")")
end

# -------------------------------------------------------------
# BEE model
# -------------------------------------------------------------
struct BeeModel <: BeeObject
    name::String
    intdict::Dict{String, BeeInt}
    booldict::Dict{String, BeeBool}
    conslist::Array{BeeConstraint, 1}
end
BeeModel(name::String) = BeeModel(name, Dict{String, BeeInt}(), Dict{String, BeeBool}(), Array{BeeConstraint,1}())

show(io::IO, m::BeeModel) = print(io, "BEE model [$(m.name)]")
show(io::IO, ::MIME"text/plain", m::BeeModel) = print(io,
"""BEE model [$(m.name)]:
* Integer variables: $(collect(values(m.intdict)))
* Boolean variables: $(collect(values(m.booldict)))
* Constraint: $(m.conslist)""")

"Check if the model has a variable called `name`"
hasvar(model::BeeModel, name::String) = haskey(model.intdict, name) || haskey(model.booldict, name)

"Check if the model has a bollean variable called `name` in `model`"
hasbool(model::BeeModel, name) = haskey(model.booldict, name)

"Check if the model has a bollean variable called `name` in `GBL_MODEL`"
hasbool(name) = hasbool(GBL_MODEL, name)

"Check if the model has a integer variable called `name` in `model`"
hasint(model::BeeModel, name) = haskey(model.intdict, name)

"Check if the model has a integer variable called `name` in `GBL_MODEL`"
hasint(name) = hasint(GBL_MODEL, name)

"Retrive an existing integer variable called `name` in `GBL_MODEL`"
getint(model, name) = model.intdict[name]

"Retrive an existing integer variable called `name` in `GBL_MODEL`"
getint(name) = getint(GBL_MODEL, name)

"Either create or retrive an existing boolean variable called `name` in `model`"
function fetchbool(model::BeeModel, name)
    if hasbool(model, name)
        model.booldict[name]
    else
        beebool(model, name)
    end
end

"Either create or retrive an existing boolean variable called `name` in `GBL_MODEL`"
fetchbool(name) = fetchbool(GBL_MODEL, name)

"Retrive an existing boolean variable called `name` in `GBL_MODEL`"
getbool(model, name) = model.booldict[name]

"Retrive an existing boolean variable called `name` in `GBL_MODEL`"
getbool(name) = getbool(GBL_MODEL, name)

"Return the number of boolean variables in `GBL_MODEL`"
boolnum() = length(GBL_MODEL.booldict)

"Return the number of integer variables in `GBL_MODEL`"
intnum() = length(GBL_MODEL.intdict)

"""
    constrain(model, cons)

Add the `cons` to `model`.  Note that unlike a variable, a constraint is not automatically added to
any model when it is created.
"""
function constrain(model::BeeModel, cons::BeeConstraint)
    push!(model.conslist, cons)
    cons
end
constrain(cons) = constrain(GBL_MODEL, cons)

macro constrain(cons)
    :(constrain($(esc(cons))))
end

"Render the global model `GBL_MODEL` to BEE syntax and print it to `stdout`."
render() = render(Base.stdout, GBL_MODEL)

"Render the global model `GBL_MODEL` to BEE syntax and print it to `io`."
render(io::IO) = render(io, GBL_MODEL)

function render(io::IO, model::BeeModel)
    for intv in values(model.intdict)
        render(io, intv)
    end
    for boolv in values(model.booldict)
        render(io, boolv)
    end
    for cons in model.conslist
        render(io, cons)
    end
    print(io, "solve satisfy\n")
end

"Delete all variables and constraints from the default model."
function reset()
    global GBL_MODEL = BeeModel("defaul model")
end

# -------------------------------------------------------------
# Call BumbleBEE directly from Julia
# -------------------------------------------------------------
struct BeeSolution
    sat::Bool
    intdict::Dict{String, Int}
    booldict::Dict{String, Bool}
end

"Call `BumbleBEE` to solve the `model` and print the output into `io`"
function solve(model::BeeModel, io::IO, deletefile=true)
    # Find where is BumbleBEE
    beepath = Sys.which("BumbleBEE")
    beedir = dirname(beepath)

    # Render solution to the file
    tempf = tempname() * ".bee"
    open(tempf, "w") do io
        render(io)
    end

    tempout = tempname() * ".sol"

    # Solve with BumbleBEE
    curdir = pwd()
    cd(beedir)
    println(tempout)
    run(pipeline(`./BumbleBEE $tempf`, stdout=tempout))
    output = readlines(tempout)
    cd(curdir)

    if deletefile
        rm(tempf)
        rm(tempout)
    end

    # filter comments
    rc = r"^%"
    ri = r"^\s*(\w+)\s*=\s*(-?+\d+)"
    rb = r"^\s*(\w+)\s*=\s*(true|false)"
    runsat = r"=====UNSATISFIABLE====="

    sat = true

    intdict = Dict{String, Int}()
    booldict = Dict{String, Bool}()

    for line in output
        println(io, line)
        if match(runsat, line) !== nothing
            sat = false
        end
        if match(rc, line) !== nothing
            continue
        elseif (m  = match(ri, line)) !== nothing
            name, val = m.captures
            intdict[name] = parse(Int, val)
        elseif (m  = match(rb, line)) !== nothing
            name, val = m.captures
            booldict[name] = parse(Bool, val)
        end
    end

    BeeSolution(sat, intdict, booldict)
end

" Solve the default model and print the solution to `stdout`."
solve() = solve(GBL_MODEL, Base.stdout)

solve(deletefile::Bool) = solve(GBL_MODEL, Base.stdout, deletefile)

show(io::IO, ::MIME"text/plain", sol::BeeSolution) = print(io,
"""
BEE solution:
* Satisfiable: $(sol.sat)
* Integer variables: $(sol.intdict)
* Boolean variables: $(sol.booldict)""")

# -------------------------------------------------------------
# BEE operator for both integer and boolean
# -------------------------------------------------------------

# Create BEE summation expression, which applies to list of symbols
struct BeeSum{T <: BeeSymbol} <: BeeExpression
    varlist::Array{BeeSymbol, 1}
end

sum(varlist::Array{T, 1})  where T <: BeeBoolean = BeeSum{T}(varlist)
sum(varlist::Array{T, 1})  where T <: BeeInteger = BeeSum{T}(varlist)

# Create BEE operator on summation.
for (VT, VF) in [(:BeeBoolean, :bool), (:BeeInteger, :int)],
    (OP, EF) in [(:(<=), :leq), (:(>=), :geq), (:(==), :eq), (:(<), :lt), (:(>), :gt), (:(!=), :neq)]
    @eval BeeEncoder begin
        # Some of these are not symmetirc. Let's don't switch order
        # $OP(lhs::ZZ, rhs::BeeSum{T} where T <: $VT)  = $OP(rhs, lhs)
        $OP(lhs::BeeSum{T} where T <: $VT, rhs::ZZ) = $(Symbol(VF, :_array_sum_, EF))(lhs.varlist, rhs)
    end
end

# -------------------------------------------------------------
# BEE operator for integers
# -------------------------------------------------------------

# Boolean operator on two integers
intBOP = [(:BeeLeq, :leq, :(<=)), (:BeeGeq, :geq, :(>=)), (:BeeEq, :eq, :(==)),
          (:BeeLt, :lt, :(<)), (:BeeGt, :gt, :(>)), (:BeeNeq, :neq, :(!=))]
# Arithmetic operator on two integers
intAOP = [(:BeePlus, :plus, :+), (:BeeTimes, :times, :*), (:BeeMax, :max, :max),
          (:BeeMin, :min, :min), (:BeeDiv, :div, :/), (:BeeMod, :mod, :mod)]

# Define function for integer $OP integer. Avoid matching `Int` $OP `Int`
for (ET, EF, OP) in [intBOP; intAOP]
    @eval BeeEncoder begin
    $OP(lhs::Int, rhs::BeeInteger) = $ET(BeeIntLiteral(lhs), rhs)
    $OP(lhs::BeeInteger, rhs::Int) = $ET(lhs, BeeIntLiteral(rhs))
    $OP(lhs::BeeInteger, rhs::BeeInteger) = $ET(lhs, rhs)
    end
end

# Create BEE boolean expression for two integers, which applies to two `BeeInteger`.
for (ET, EF, OP) in intBOP
    @eval BeeEncoder begin
    struct $ET <: BeeExpression
        lhs::BeeInteger
        rhs::BeeInteger
    end

    # `lhs` is true `iff` rhs is true
    (==)(lhs::BeeBoolean, rhs::$ET) = rhs == lhs
    function (==)(lhs::$ET, rhs::BeeBoolean)
        $(Symbol(:int_, EF, :_reif))(lhs.lhs, lhs.rhs, rhs)
    end

    # `lhs` is true `iff` rhs is true
    (==)(lhs::Bool, rhs::$ET) = rhs == lhs
    function (==)(lhs::$ET, rhs::Bool)
        if rhs
            $(Symbol(:int_, EF))(lhs.lhs, lhs.rhs)
        else
            $(Symbol(:int_, EF, :_reif))(lhs.lhs, lhs.rhs, rhs)
        end
    end
    end
end

# Create BEE arithmetic expression for two integers, which applies to two `BeeInteger`.
for (ET, EF, OP) in  intAOP
    @eval BeeEncoder begin
    struct $ET <: BeeExpression
        lhs::BeeInteger
        rhs::BeeInteger
    end

    # `lhs` == `rhs` is true
    (==)(lhs::ZZ, rhs::$ET) = rhs == lhs
    function (==)(lhs::$ET, rhs::ZZ)
        $(Symbol(:int_, EF))(lhs.lhs, lhs.rhs, rhs)
    end
    end
end

# Create BEE allDiff operations on one integer array
intarrayOP = [(:BeeArrayAllDiff, :allDiff, :alldiff, :_reif)]
for (ET, EF, OP, TAIL) in intarrayOP
    @eval BeeEncoder begin
    struct $ET <: BeeExpression
        varlist::Array{BeeInteger, 1}
    end

    # No need to check type here
    $OP(varlist::Array{T, 1} where T <: ZZ) = $ET(varlist)

    (==)(lhs::BeeInteger, rhs::$ET) = rhs == lhs
    function (==)(lhs::$ET, rhs::BeeInteger)
        $(Symbol(:int_array_, EF, TAIL))(lhs.varlist, rhs)
    end

    (==)(lhs::Bool, rhs::$ET) = rhs == lhs
    function (==)(lhs::$ET, rhs::Bool)
        if rhs
            $(Symbol(:int_array_, EF))(lhs.varlist)
        else
            $(Symbol(:int_array_, EF, TAIL))(lhs.varlist, rhs)
        end
    end
    end
end

# Create BEE operations on one integer array
intarrayOP = [(:BeeArrayPlus, :plus, :plus),
              (:BeeArrayTimes, :times, :times),
              (:BeeArrayMax, :max, :max),
              (:BeeArrayMin, :min, :min)]
for (ET, EF, OP) in intarrayOP
    @eval BeeEncoder begin
    struct $ET <: BeeExpression
        varlist::Array{BeeInteger, 1}
    end

    # No need to check type here
    $OP(varlist::Array{T, 1} where T <: ZZ) = $ET(varlist)

    (==)(lhs::BeeInteger, rhs::$ET) = rhs == lhs
    function (==)(lhs::$ET, rhs::BeeInteger)
        $(Symbol(:int_array_, EF))(lhs.varlist, rhs)
    end
    end
end

# -------------------------------------------------------------
# BEE operator for boolean
# -------------------------------------------------------------

# hard code this bit since it's just one function
(==)(lhs::Bool, rhs::BeeBoolean) = bool_eq(lhs, rhs)
(==)(lhs::BeeBoolean, rhs::Bool) = bool_eq(lhs, rhs)
(==)(lhs::BeeBoolean, rhs::BeeBoolean) = bool_eq(lhs, rhs)

# Logic Expressions on two boolean.
boolOP = [(:BeeAnd, :and, :and), (:BeeOr, :or, :or), (:BeeXor, :xor, :xor), (:BeeIff, :iff, :iff)]
for (ET, EF, OP) in boolOP
    @eval BeeEncoder begin
    struct $ET <: BeeExpression
        lhs::BeeSymbol
        rhs::BeeSymbol
    end

    # No need to check type here
    $OP(lhs, rhs) = $ET(lhs, rhs)

    (==)(lhs::BeeBoolean, rhs::$ET) = rhs == lhs
    function (==)(lhs::$ET, rhs::BeeBoolean)
        $(Symbol(:bool_, EF, :_reif))(lhs.lhs, lhs.rhs, rhs)
    end

    (==)(lhs::Bool, rhs::$ET) = rhs == lhs
    function (==)(lhs::$ET, rhs::Bool)
        if rhs
            $(Symbol(:bool_, EF))(lhs.lhs, lhs.rhs)
        else
            $(Symbol(:bool_, EF, :_reif))(lhs.lhs, lhs.rhs, rhs)
        end
    end
    end
end

# Create BEE operators on boolean arrays. Note that
boolarrayOP = [(:BeeArrayAnd, :and, :and), (:BeeArrayOr, :or, :or), (:BeeArrayXor, :xor, :xor), (:BeeArrayIff, :iff, :iff)]
for (ET, EF, OP) in boolarrayOP
    @eval BeeEncoder begin
    struct $ET <: BeeExpression
        varlist::Array{BeeBoolean, 1}
    end

    # No need to check type here
    $OP(varlist) = $ET(varlist)

    (==)(lhs::BeeBoolean, rhs::$ET) = rhs == lhs
    function (==)(lhs::$ET, rhs::BeeBoolean)
        $(Symbol(:bool_array_, EF, :_reif))(lhs.varlist, rhs)
    end

    (==)(lhs::Bool, rhs::$ET) = rhs == lhs
    function (==)(lhs::$ET, rhs::Bool)
        if rhs
            $(Symbol(:bool_array_, EF))(lhs.varlist)
        else
            $(Symbol(:bool_array_, EF, :_reif))(lhs.varlist, rhs)
        end
    end
    end
end

# Create BEE operators on two boolean arrays
bool2arrayOP = [(:BeeArrayEq, :eq, :(==)), (:BeeArrayNeq, :neq, :(!=)),
                     (:BeeLex, :lex, :(<=)), (:BeelexLt, :lexLt, :(<))]
for (ET, EF, OP) in bool2arrayOP
    @eval BeeEncoder begin
    struct $ET <: BeeExpression
        lhs::Array{BeeBoolean, 1}
        rhs::Array{BeeBoolean, 1}
    end

    $OP(lhs::Array{T, 1} where T <: BB, rhs::Array{T, 1} where T <: BB) = $ET(lhs, rhs)

    (==)(lhs::BeeBoolean, rhs::$ET) = rhs == lhs
    function (==)(lhs::$ET, rhs::BeeBoolean)
        $(Symbol(:bool_arrays_, EF, :_reif))(lhs.lhs, lhs.rhs, rhs)
    end

    (==)(lhs::Bool, rhs::$ET) = rhs == lhs
    function (==)(lhs::$ET, rhs::Bool)
        if rhs
            $(Symbol(:bool_arrays_, EF))(lhs.lhs, lhs.rhs)
        else
            $(Symbol(:bool_arrays_, EF, :_reif))(lhs.lhs, lhs.rhs, rhs)
        end
    end
    end
end

# -------------------------------------------------------------
# BEE functions
# -------------------------------------------------------------

# Adding constraint function. Do not do any check. Let BEE to report errors.
for F in (:int_order2bool_array,
          :bool2int,
          :bool_eq,
          :bool_array_eq_reif,
          :bool_array_or, :bool_array_and, :bool_array_xor, :bool_array_iff,
          :bool_array_or_reif, :bool_array_and_reif, :bool_array_xor_reif, :bool_array_iff_reif,
          :bool_or_reif, :bool_and_reif, :bool_xor_reif, :bool_iff_reif,
          :bool_ite,
          :bool_ite_reif,
          :int_leq, :int_geq, :int_eq, :int_lt, :int_gt, :int_neq,
          :int_leq_reif, :int_geq_reif, :int_eq_reif, :int_lt_reif, :int_gt_reif, :int_neq_reif,
          :int_array_allDiff,
          :int_array_allDiff_reif,
          :int_array_allDiffCond,
          :int_abs,
          :int_plus, :int_times, :int_div, :int_mod, :int_max, :int_min,
          :int_array_plus, :int_array_times, :int_array_max, :int_array_min,
          :bool_array_sum_leq, :bool_array_sum_geq, :bool_array_sum_eq, :bool_array_sum_lt, :bool_array_sum_gt,
          :bool_array_pb_leq, :bool_array_pb_geq, :bool_array_pb_eq, :bool_array_pb_lt, :bool_array_pb_gt,
          :int_array_sum_leq, :int_array_sum_geq, :int_array_sum_eq, :int_array_sum_lt, :int_array_sum_gt,
          :int_array_lin_leq, :int_array_lin_geq, :int_array_lin_eq, :int_array_lin_lt, :int_array_lin_gt,
          :int_array_sumCond_leq, :int_array_sumCond_geq, :int_array_sumCond_eq, :int_array_sumCond_lt, :int_array_sumCond_gt,
          :bool_array_sum_modK,
          :bool_array_sum_divK,
          :bool_array_sum_divModK,
          :int_array_sum_modK,
          :int_array_sum_divK,
          :int_array_sum_divModK,
          :bool_arrays_eq, :bool_arrays_neq,
          :bool_arrays_eq_reif, :bool_arrays_neq_reif,
          :bool_arrays_lex,
          :bool_arrays_lexLt,
          :bool_arrays_lex_reif,
          :bool_arrays_lexLt_reif,
          :int_arrays_eq, :int_arrays_neq,
          :int_arrays_eq_reif, :int_arrays_neq_reif,
          :int_arrays_lex,
          :int_arrays_lexLt,
          :int_arrays_lex_implied,
          :int_arrays_lexLt_implied,
          :int_arrays_lex_reif,
          :int_arrays_lexLt_reif)
    SF = String(F)
    @eval BeeEncoder $F(var...) = BeeConstraint($SF, [var...])
end

# -------------------------------------------------------------
# Initialize module
# -------------------------------------------------------------

reset()

# -------------------------------------------------------------
# For testing
# -------------------------------------------------------------
function capture_render(c)
    @capture_out render(c)
end

function capture_render()
    @capture_out render()
end

end

max_stars_repo_name: newptcai/BeeEncoder

max_issues_repo_name: newptcai/BeeEncoder

max_forks_repo_name: newptcai/BeeEncoder


text:
instructions = []
for line in eachline(STDIN, chomp = true)
  push!(instructions, parse(Int, line))
end
len = length(instructions)

index = 1
steps = 0
while index > 0 && index <= len
  jump = instructions[index]
  instructions[index] += (jump >= 3) ? -1 : 1
  index += jump
  steps += 1
end

println("Step count for part 2: $steps")


text:
function register!(actr::AbstractACTR, fun, when::Now, args...; id="", type="", description="", kwargs...)
    scheduler = actr.scheduler
    register!(scheduler, fun, scheduler.time, args...; id, type, description, kwargs...)
end

function register!(actr::AbstractACTR, fun, when::At, t, args...; id="", type="", description="", kwargs...)
    scheduler = actr.scheduler
    register!(scheduler, fun, t, args...; id, type, description, kwargs...)
end

function register!(actr::AbstractACTR, fun, when::After, t, args...; id="", type="", description="", kwargs...)
    scheduler = actr.scheduler
    register!(scheduler, fun, scheduler.time + t, args...; id, type, description, kwargs...)
end

function register!(actr::AbstractACTR, fun, when::Every, t, args...; id="", type="", description="", kwargs...)
    scheduler = actr.scheduler
    function f(args...; kwargs...)
        fun1 = ()->fun(args...; kwargs...)
        fun1()
        register!(scheduler, fun, every, t, args...; id, type, description, kwargs...)
    end
    register!(scheduler, f, after, t, args...; id, type, description, kwargs...)
end

function register!(task::AbstractTask, fun, when::Now, args...; id="", type="", description="", kwargs...)
    scheduler = task.scheduler
    register!(scheduler, fun, scheduler.time, args...; id, type, description, kwargs...)
end

function register!(task::AbstractTask, fun, when::At, t, args...; id="", type="", description="", kwargs...)
    scheduler = task.scheduler
    register!(scheduler, fun, t, args...; id, type, description, kwargs...)
end

function register!(task::AbstractTask, fun, when::After, t, args...; id="", type="", description="", kwargs...)
    scheduler = task.scheduler
    register!(scheduler, fun, scheduler.time + t, args...; id, type, description, kwargs...)
end

function register!(task::AbstractTask, fun, when::Every, t, args...; id="", type="", description="", kwargs...)
    scheduler = task.scheduler
    function f(args...; kwargs...)
        fun1 = ()->fun(args...; kwargs...)
        fun1()
        register!(scheduler, fun, every, t, args...; id, type, description, kwargs...)
    end
    register!(scheduler, f, after, t, args...; id, type, description, kwargs...)
end

function import_gui()
    path = pathof(ACTRSimulators) |> dirname |> x->joinpath(x, "")
    include(path * "GUI.jl")
end

function press_key!(task, actr, key)
    @error "A method must be defined for press_key! in the form of
    press_key!(task, actr, key)"
end

function start!(task, actr)
    @error "A method must be defined for start! in the form of
    start!(task, actr)"
end

"""
    why_not(actr, rule)

Prints each condition of a production rule and whether it matches the state of the
architecture.

# Arguments

- `actr`: an ACT-R model object
- `rule`: a production rule object
"""
function why_not(actr, rule)
    str = rule.name * "\n"
    for c in rule.conditions
         str *= @code_string c(actr)
         str *= string(" ", c(actr))
         str *= "\n"
    end
    println(str)
    return nothing
end

"""
    why_not(actr)

Prints each condition of all production rules and whether the condition matche the state of the
architecture.

# Arguments

- `actr`: an ACT-R model object
"""
function why_not(actr)
    return why_not.(actr, get_rules(actr))
end


text:
"""
    datasetpath(name)

Return the folder that dataset `name` is stored.

If it hasn't been downloaded yet, you will be asked if you want to
download it. See [`Datasets.DATASETS`](#) for a list of available datasets.
"""
function datasetpath(name)
    i = findfirst(DATASETS .== name)
    isnothing(i) && error("Dataset $name does not exist. Check `DATASETS` for available datasets.")
    config = DATASETCONFIGS[i]
    datadeppath = @datadep_str "fastai-$name"
    return Path(joinpath(datadeppath, config.subpath))
end

function loadfolderdata(
        dir;
        pattern="**",
        splitfn = nothing,
        filterfn = nothing,
        loadfn = nothing)
    data = FileDataset(dir, pattern)
    if filterfn !== nothing
        data = filterobs(filterfn, data)
    end
    if splitfn !== nothing
        data = groupobs(splitfn, data)
    end
    if loadfn !== nothing
        if splitfn === nothing
            data = mapobs(loadfn, data)
        else
            data = Dict(zip(keys(data), map(d -> mapobs(loadfn, d), values(data))))
        end
    end
    return data
end

parentname(f) = f |> pathparent |> pathname
grandparentname(f) = f |> pathparent |> pathparent |> pathname
matches(re::Regex) = f -> matches(re, f)
matches(re::Regex, f) = !isnothing(match(re, f))
const RE_IMAGEFILE = r".*\.(gif|jpe?g|tiff?|png|webp|bmp)$"i
isimagefile(f) = matches(RE_IMAGEFILE, f)

maskfromimage(a::AbstractArray{<:Gray{T}}, classes) where T = maskfromimage(reinterpret(T, a), classes)
maskfromimage(a::AbstractArray{<:Normed{T}}, classes) where T = maskfromimage(reinterpret(T, a), classes)
function maskfromimage(a::AbstractArray{I}, classes) where {I<:Integer}
    a .+= one(I)
    return IndirectArray(a, classes)
end

"""
    loadmask(file, classes)

Load a segmentation mask from an image file. Returns an efficiently stored
array of type `eltype(classes)`.

"""
loadmask(file, classes) = maskfromimage(loadfile(file), classes)


text:
#! /usr/bin/env julia

using Docile, # using rewrite branch as of 2015-11-26
      GeometryTypes

makedocs()


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/functions/nonlinear/functions.jl
# Defines variable info and function types.
#############################################################################

"""
$(TYPEDEF)

An abstract super-type used for representing constraints built by EAGO's backend.
"""
abstract type AbstractEAGOConstraint end

"""

Computes the lower interval bound for `AbstractEAGOConstraint` representing an
inequality constraint.
"""
function lower_interval_bound end

"""

Computes a tuple representing the lower and upper interval bounds for a
`AbstractEAGOConstraint` representing an equality constraint.
"""
function interval_bound end

"""

Eliminate fixed variables by rearrangment or restructuring of `AbstractEAGOConstraint`.
"""
function eliminate_fixed_variables! end

include("affine.jl")
include("quadratic.jl")
include("cone.jl")
include("nonlinear.jl")


text:
# A comparison between ℓ_0 and ℓ_1 trend filtering

using Plots
using Convex
using Mosek
using EllZeroTrendFiltering

g = snp500_data()[1:300]

N = length(g)

## ℓ_1 trend filtering, see "ℓ_1 Trend Filtering", Kim et al. (2009)
H =  spdiagm((ones(N-2), -2*ones(N-2), ones(N-2)), (0,1,2))

x = Variable(N)

λ = 2.0
problem = minimize(sumsquares(x - g) + λ * norm(H*x, 1))

solve!(problem, MosekSolver())
println(problem.status)

x_l1 = evaluate(x)[:]
f_l1 = sum((x_l1 - g).^2)

## ℓ_0 trend filtering
ζ = 0.02
I_l0, Y_l0, f_l0 = fit_pwl_regularized(g, ζ)

plot(g, label="SP500", size=(1000,550));
plot!(x_l1, lw=2.5, label="l1 reg. lambda = $λ");
plot!(I_l0, Y_l0, lw=2.5, label="l0 reg. zeta = $ζ");
gui()

plot(H*x_sol);
plot!(I_vec[m][2:end], 0.1*diff(Y_vec[m]) ./ diff(I_vec[m]), m=:o, lw=0)


text:
using DebuggerFramework
using Cxx
using Iterators
reload("DebuggerFramework")
module LLVMIRDebugger
  using DebuggerFramework
  using Cxx
  using Iterators

  Cxx.addHeaderDir("/home/keno/julia/deps/srccache/llvm-3.9.1/")
  cxx"""
  #include "llvm/IR/Function.h"
  #include "llvm/IR/Module.h"
  #include "llvm/IR/BasicBlock.h"
  #include "llvm/IRReader/IRReader.h"
  #include "llvm/Support/MemoryBuffer.h"
  #include "llvm/IR/ModuleSlotTracker.h"

  // This is a private LLVM header at the moment
  #include "lib/ExecutionEngine/Interpreter/Interpreter.h"

  extern llvm::LLVMContext jl_LLVMContext;
  """

  # The interpreter
  struct ExecutionContext
      M::pcpp"llvm::Module"
      Interp::pcpp"llvm::Interpreter"
      ExecutionContext(M::pcpp"llvm::Module") = new(M,
          icxx"new llvm::Interpreter(std::unique_ptr<llvm::Module>{$M});")
  end
  CurBB(ctx::ExecutionContext) = icxx"$(ctx.Interp)->ECStack.back().CurBB;"
  CurInst(ctx::ExecutionContext) = icxx"$(ctx.Interp)->ECStack.back().CurInst;"
  CurFunction(ctx::ExecutionContext) = icxx"$(ctx.Interp)->ECStack.back().CurFunction;"

  function Base.convert(::Type{vcpp"llvm::MemoryBufferRef"}, data::String)
      icxx"""
        llvm::MemoryBufferRef{StringRef{$(pointer(data)),(size_t)$(sizeof(data))},
                        "converted"};
      """
  end

  function Base.parse(::Type{pcpp"llvm::Module"}, data)
      mod = icxx"""
      llvm::SMDiagnostic Err;
      auto mod =
        llvm::parseIR($(convert(vcpp"llvm::MemoryBufferRef", data)),
                      Err, jl_LLVMContext);
      if (!mod) {
        Err.print("llvmir.jl", llvm::errs());
      }
      mod.release();
      """
      if mod == C_NULL
          error("Parsing failed")
      end
      mod
  end

  function Base.print(io::IO, Val::pcpp"llvm::Value")
      print(io, unsafe_string(icxx"""
          std::string str;
          llvm::raw_string_ostream OS(str);
          $Val->print(OS, true);
          OS.str();
          str;
      """))
  end

  Base.iterate(BB::pcpp"llvm::BasicBlock", i=icxx"$BB->begin();") =
    icxx"$i == $BB->end();" ? nothing : (icxx"*$i;", icxx"++$i;")
  prev(BB::pcpp"llvm::BasicBlock", i) = (icxx"*$i;", icxx"--$i;")

  function step_one!(state::ExecutionContext)
      # Don't use the acessor, because we also want to increment
      icxx"$(state.Interp)->visit(*$(state.Interp)->ECStack.back().CurInst++);"
  end

  # Debugger stuff
  function enter(F::pcpp"llvm::Function")
      ec = ExecutionContext(icxx"$F->getParent();")
      icxx"$(ec.Interp)->callFunction($F, llvm::ArrayRef<llvm::GenericValue>());"
      ec
  end

  # Debugger support
  struct LLVMIRFrame <: DebuggerFramework.StackFrame
      state::ExecutionContext
      idx::Int
  end
  CurFunction(frame::LLVMIRFrame) = icxx"$(frame.state.Interp)->ECStack[$(frame.idx)].CurFunction;"

  function DebuggerFramework.debug(F::pcpp"llvm::Function", args...)
      state = enter(F)
      stack = [LLVMIRFrame(state, 0)]
      DebuggerFramework.RunDebugger(stack, args...)
  end

  function getName(V::Union{pcpp"llvm::BasicBlock", pcpp"llvm::Function", pcpp"llvm::Value"})
      unsafe_string(icxx"$(V)->getName();")
  end

  function DebuggerFramework.print_status_synthtic(io::IO, state, frame::LLVMIRFrame, lines_before, total_lines)
      es = frame.state
      cur_it = icxx"auto it = $(CurInst(es)); it;" # Fake copy constructor
      show_bb_name = false
      for i = 1:lines_before
          if icxx"&*$(cur_it) == &$(CurBB(es))->front();"
              show_bb_name = true
              break
          end
          cur_it = prev(CurBB(es), cur_it)[2]
      end
      active_line = cur_line = 0
      if show_bb_name
          println(io, getName(CurBB(es)),":")
          cur_line += 1
          total_lines -= 1
      end
      for inst in Base.Iterators.take(Rest(CurBB(es), cur_it), 5)
          cur_line += 1
          if icxx"&$inst == &*$(CurInst(es));"
              active_line = cur_line
          end
          println(io, icxx"(llvm::Value*)&$inst;")
      end
      active_line
  end

  function DebuggerFramework.locdesc(frame::LLVMIRFrame)
      F = CurFunction(frame)
      BB = icxx"$(frame.state.Interp)->ECStack[$(frame.idx)].CurBB;"
      "function $(getName(F)), BB $(getName(BB))"
  end

  struct AggregateValue
      values::Vector{Any}
  end
  Base.print(io::IO, val::AggregateValue) = print(io, val.values)

  function fromGenericValue(typ::pcpp"llvm::Type", val::vcpp"llvm::GenericValue")
      if icxx"$typ->isVoidTy();"
          return nothing
      elseif icxx"$typ->isHalfTy() || $typ->isX86_FP80Ty() ||
                  $typ->isFP128Ty() || $typ->isPPC_FP128Ty();"
          error("Unsupported floating point type")
      elseif icxx"$typ->isFloatTy();"
          return icxx"$val.FloatVal;"
      elseif icxx"$typ->isDoubleTy();"
          return icxx"$val.DoubleVal;"
      elseif icxx"$typ->isPointerTy();"
          return icxx"$val.PointerVal;"
      elseif icxx"$typ->isAggregateType();"
          return AggregateValue(Any[fromGenericValue(x...) for x in
            zip(icxx"$typ->subtypes();", icxx"$val.AggregateVal;")])
      elseif icxx"$typ->isIntegerTy();"
          # Good enough for now
          return parse(BigInt, unsafe_string(icxx"$val.IntVal.toString(10, false);"))
      else
          error("Unsupported type")
      end
  end

  function DebuggerFramework.print_locals(io::IO, frame::LLVMIRFrame)
      MST = icxx"new llvm::ModuleSlotTracker{$(frame.state.M), true};"
      icxx"$MST->incorporateFunction(*$(CurFunction(frame)));"
      for (CompileTimeVal, RunTimeVal) in icxx"$(frame.state.Interp)->ECStack[$(frame.idx)].Values;"
          name = string("%",icxx"$CompileTimeVal->hasName();" ?
            getName(CompileTimeVal) :
            string(icxx"$MST->getLocalSlot($CompileTimeVal);"))
          DebuggerFramework.print_var(io, name, Nullable(
              fromGenericValue(icxx"$CompileTimeVal->getType();", RunTimeVal)
            ), nothing)
      end
      icxx"delete $MST;"
  end

  function update_stack!(ds, state)
      stacksize = Int(icxx"$(state.Interp)->ECStack.size();"-1)
      ds.stack = [LLVMIRFrame(state, i) for i = stacksize:-1:0]
      if ds.level > length(ds.stack)
          ds.level = length(ds.stack)
      end
  end

  function DebuggerFramework.execute_command(state, frame::LLVMIRFrame, cmd::Val{:si}, command)
      step_one!(frame.state)
      update_stack!(state, frame.state)
      return true
  end
end

ir = """
@.str = private unnamed_addr constant [13 x i8] c"hello world\\0A\\00"
declare i32 @puts(i8* nocapture) nounwind
define void @julia_foo(i8*, i8*) {
top:
  br label %if

if:                                               ; preds = %top, %pass
  %"#temp#.02" = phi i64 [ 1, %top ], [ %4, %if ]
  %2 = and i64 %"#temp#.02", 255
  %3 = icmp eq i64 %2, %"#temp#.02"
  %4 = add nuw nsw i64 %"#temp#.02", 1
  %5 = trunc i64 %"#temp#.02" to i8
  %6 = add nsw i64 %"#temp#.02", -1
  %7 = getelementptr i8, i8* %1, i64 %6
  %8 = load i8, i8* %7, align 1
  %9 = zext i8 %8 to i64
  %10 = add nsw i64 %9, -1
  %11 = getelementptr i8, i8* %0, i64 %10
  store i8 %5, i8* %11, align 1
  %12 = icmp eq i64 %4, 21
  br i1 %12, label %done, label %if

done:
  ret void
}
define i32 @main() {   ; i32()*
top:
    %temp = alloca [13 x i8]
    %casttemp = getelementptr [13 x i8], [13 x i8]* %temp, i64 0, i64 0
    %cast210 = getelementptr [13 x i8], [13 x i8]* @.str, i64 0, i64 0
    call void @julia_foo(i8* %casttemp, i8* %cast210)
    call i32 @puts(i8* %cast210)
    ret i32 0
}
"""
M = parse(pcpp"llvm::Module", ir)
F = icxx"""$M->getFunction("main");"""


text:
using Documenter, ComputerVision

makedocs(;
    modules=[ComputerVision],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/zygmuntszpak/ComputerVision.jl/blob/{commit}{path}#L{line}",
    sitename="ComputerVision.jl",
    authors="Dr. Zygmunt L. Szpak",
    assets=String[],
)

deploydocs(;
    repo="github.com/zygmuntszpak/ComputerVision.jl",
)


text:
export timeseries, segment

"""
    timeseries(notes::Notes, property::Symbol, f, grid; kwargs...) -> tvec, ts
Produce a timeseries of the `property` of the given notes, by
first quantizing on the given `grid` (to avoid actual quantization
use the grid `0:1//notes.tpq:1`). Return the time vector `tvec` in ticks
and the produced timeseries `ts`.

After quantization, it is often the case that many notes are in the same bin of
the grid. The **function** `f` denotes which value of the vector of the
`property` of the notes to keep. Typical values are `minimum, maximum, mean`,
etc. Notice that bins without any note in them obtain the value of the keyword `missingval`,
which be default is just `missing`, regardless of the function `f` or the `property`.

If the `property` is `:velocity`, `:pitch`, or `:duration` the function
behaves exactly as described. The `property` can also be `:position`.
In this case, the timeseries `ts` contain the timing deviations of the notes
with respect to the `tvec` vector
(these numbers are known as *microtiming deviations* in the literature).

If given keyword 'segmented = true', the notes are segmented according to the grid
in order to respect the information of their duration, see [`segment`](@ref).
Otherwise the notes are treated as point events with no duration
(it makes no sense to choose `:duration` with `segmented`).
"""
function timeseries(notes, property, f, grid; segmented = false, missingval = missing)
    isgrid(grid)
    if segmented == true
        notes = segment(notes, grid)
    end
    if !issorted(notes, by = x -> x.position)
        error("notes must be sorted by position!")
    elseif !isnothing(property) &&
        property ∉ (:velocity, :pitch, :duration, :position, :channel)
        error("Unknown property!")
    end
    ts, tvec, quantizedpos = _init_timeseries_vectors(notes, grid, missingval)
    i = previdx = 1; L = length(quantizedpos); M = length(tvec)
    while i ≤ L
        # find entries of same grid bin
        j = 1
        while j ≤ L - i && quantizedpos[i+j] == quantizedpos[i]
            j+=1
        end
        # where to add the value in the timeseries:
        idx = findfirst(x -> tvec[x] == quantizedpos[i], previdx:M) + previdx-1
        previdx = idx
        ts[idx] = produce_note_value(notes, property, f, i, j)
        if property == :position # here we want timing *deviations*
            ts[idx] -= tvec[idx]
        end
        i += j
    end
    return tvec, ts
end

"""
    segment(notes, grid) → segmented_notes
Quantize the positions and durations of `notes` and then segment them (i.e. cut them into
pieces) according to the duration of a grid unit. This function only works with
`AbstractRange` grids, i.e. equi-spaced grids like `0:1//3:1`.
"""
function segment(notes::Notes{N}, grid::AbstractRange) where {N}
    isgrid(grid)
    segment_duration = Int(step(grid)*notes.tpq)
    qnotes = quantize(notes, grid)
    segmented_notes = Notes(Vector{N}(), notes.tpq)
    for qn in qnotes
        for segment_idx in 0:round(Int, qn.duration/segment_duration)-1
            push!(segmented_notes, N(qn.pitch, qn.velocity,
                  qn.position + segment_idx*segment_duration, segment_duration))
        end
    end
    return timesort!(segmented_notes)
end

"""
    timeseries(notes::Notes, f, grid) -> tvec, ts
If `property` is not given, then `f` should take as input a `Notes` instance
and output a numeric value. This is useful for example in cases where one
would want the timeseries of the velocities of the notes of the
highest pitch.
"""
timeseries(notes, f, grid; segmented = false) = timeseries(notes, nothing, f, grid; segmented = false)

function produce_note_value(notes, property::Symbol, f, i, j)
    if j > 1
        val = Float64(f(getfield(notes[k], property) for k in i:i+j-1))
    else
        val = Float64(getfield(notes[i], property))
    end
    return val
end

function produce_note_value(notes, property::Nothing, f, i, j)
    v = notes[i:i+j-1]
    val = Float64(f(v))
end

function _init_timeseries_vectors(notes, grid, missingval)
    tpq = notes.tpq
    qnotes = quantize(notes, grid)
    quantizedpos = positions(qnotes)
    bins = round.(Int, tpq .* grid)[1:end-1]
    # tvec limits
    firstp = Int(qnotes[1].position)
    firstbin = findfirst(x -> x == mod(firstp, tpq), bins)

    # Create vectors
    tvec = (firstp÷tpq)*tpq .+ bins[firstbin:end]
    c = firstp÷tpq + 1
    while tvec[end] < notes[end].position
        append!(tvec, bins .+ c*tpq)
        c += 1
    end
    ts = Vector{Union{Float64, typeof(missingval)}}(undef, length(tvec))
    fill!(ts, missingval)
    return ts, tvec, quantizedpos
end


text:
using NDReducibles
using Referenceables

using Transducers
using Transducers: maybe_usesimd, BottomRF, SideEffect

_ndrmul_reducible(C, A, B) =
    ndreducible(
        referenceable(C) => (:i, :j),
        A => (:i, :k),
        B => (:k, :j),
    )

function _ndrmul_rf()
    simd = Val(:ivdep)

    rf = SideEffect() do (c, a, b)
        c[] = muladd(a, b, c[])
        return
    end

    return maybe_usesimd(BottomRF{Any}(rf), simd)
end

function ndrmul!(C, A, B)
    fill!(C, 0)

    CAB = _ndrmul_reducible(C, A, B)
    rf = _ndrmul_rf()
    transduce(rf, nothing, CAB)

    return C
end

function manmul!(C, A, B)
    fill!(C, 0)
    for j in 1:size(C, 2), k in 1:size(A, 2)
        b = @inbounds B[k, j]
        aofs = LinearIndices(A)[1, k] - 1
        cofs = LinearIndices(C)[1, j] - 1
        @simd for i in 1:size(A, 1)
            c = @inbounds C[i + cofs]
            a = @inbounds A[i + aofs]
            @inbounds C[i + cofs] = muladd(a, b, c)
        end
    end
    C
end


text:
import QuantumControlBase
using QuantumControlBase:
    getcontrols, getcontrolderivs, discretize_on_midpoints, evalcontrols
using QuantumControlBase: GradVector, TimeDependentGradGenerator
using QuantumControlBase.QuantumPropagators: init_storage, initpropwrk
using ConcreteStructs

# GRAPE workspace (for internal use)
@concrete terse struct GrapeWrk

    # a copy of the objectives
    objectives

    # the adjoint objectives, containing the adjoint generators for the
    # backward propagation
    adjoint_objectives

    # The kwargs from the control problem
    kwargs

    # Tuple of the original controls (probably functions)
    controls

    # TODO: pulsevals0 and pulsevals1
    pulsevals::Vector{Float64}

    gradient::Vector{Float64}  # gradient for guess in iterations

    searchdirection::Vector{Float64}  # search-direction for guess in iterations

    fg_count::Vector{Int64}

    # Result object
    result

    #################################
    # scratch objects, per objective:

    # backward-propagated states
    # note: storage for fw-propagated states is in result.states
    bw_states

    # forward-propagated states (functional evaluation only)
    fw_states

    # backward-propagated grad-vectors
    bw_grad_states

    # gradients ∂τₖ/ϵₗ(tₙ)
    tau_grads::Vector{Matrix{ComplexF64}}

    # dynamical generator (normal propagation) at a particular point in time
    G

    # dynamical generator for grad-bw-propagation, time-dependent
    TDgradG

    # dynamical generator for grad-propagation at a particular point in time
    gradG

    control_derivs::Vector{Vector{Union{Function,Nothing}}}

    vals_dict

    fw_storage # backward storage array (per objective)

    fw_prop_wrk # for normal forward propagation

    grad_prop_wrk # for gradient forward propagation

    use_threads::Bool

end

function GrapeWrk(problem::QuantumControlBase.ControlProblem; verbose=false)
    use_threads = get(problem.kwargs, :use_threads, false)
    objectives = [obj for obj in problem.objectives]
    adjoint_objectives = [adjoint(obj) for obj in problem.objectives]
    controls = getcontrols(objectives)
    control_derivs = [getcontrolderivs(obj.generator, controls) for obj in objectives]
    tlist = problem.tlist
    # interleave the pulse values as [ϵ₁(t̃₁), ϵ₂(t̃₁), ..., ϵ₁(t̃₂), ϵ₂(t̃₂), ...]
    # to allow access as reshape(pulsevals0, L :)[l, n] where l is the control
    # index and n is the time index
    pulsevals = convert(
        Vector{Float64},
        reshape(
            transpose(
                hcat([discretize_on_midpoints(control, tlist) for control in controls]...)
            ),
            :
        )
    )
    kwargs = Dict(problem.kwargs)
    pulse_options = problem.pulse_options
    fg_count = zeros(Int64, 2)
    if haskey(kwargs, :continue_from)
        @info "Continuing previous optimization"
        result = kwargs[:continue_from]
        if !(result isa GrapeResult)
            # account for continuing from a different optimization method
            result = convert(GrapeResult, result)
        end
        result.iter_stop = get(problem.kwargs, :iter_stop, 5000)
        result.converged = false
        result.start_local_time = now()
        result.message = "in progress"
        pulsevals = convert(
            Vector{Float64},
            reshape(
                transpose(
                    hcat(
                        [
                            discretize_on_midpoints(control, wrk.result.tlist) for
                            control in result.optimized_controls
                        ]...
                    )
                ),
                :
            )
        )
    else
        result = GrapeResult(problem)
    end
    gradient = zeros(length(pulsevals))
    searchdirection = zeros(length(pulsevals))
    bw_states = [similar(obj.initial_state) for obj in objectives]
    fw_states = [similar(obj.initial_state) for obj in objectives]
    bw_grad_states = [GradVector(obj.initial_state, length(controls)) for obj in objectives]
    dummy_vals = IdDict(control => 1.0 for (i, control) in enumerate(controls))
    G = [evalcontrols(obj.generator, dummy_vals) for obj in objectives]
    vals_dict = [copy(dummy_vals) for _ in objectives]
    fw_storage = [init_storage(obj.initial_state, tlist) for obj in objectives]
    fw_prop_method = [
        Val(
            QuantumControlBase.get_objective_prop_method(
                obj,
                :fw_prop_method,
                :prop_method;
                problem.kwargs...
            )
        ) for obj in objectives
    ]
    bw_prop_method = [
        Val(
            QuantumControlBase.get_objective_prop_method(
                obj,
                :bw_prop_method,
                :prop_method;
                problem.kwargs...
            )
        ) for obj in objectives
    ]
    grad_prop_method = [
        Val(
            QuantumControlBase.get_objective_prop_method(
                obj,
                :grad_prop_method,
                :bw_prop_method,
                :prop_method;
                problem.kwargs...
            )
        ) for obj in objectives
    ]

    fw_prop_wrk = [
        QuantumControlBase.initobjpropwrk(
            obj,
            tlist,
            fw_prop_method[k];
            initial_state=obj.initial_state,
            info_msg="Initializing fw-prop of objective $k",
            verbose=verbose,
            kwargs...
        ) for (k, obj) in enumerate(objectives)
    ]
    TDgradG = [TimeDependentGradGenerator(obj.generator) for obj in adjoint_objectives]
    gradG = [evalcontrols(G̃_of_t, dummy_vals) for G̃_of_t ∈ TDgradG]
    grad_prop_wrk = [
        _init_objgradpropwrk(
            Ψ̃,
            tlist,
            grad_prop_method[k],
            gradG[k];
            verbose=verbose,
            info_msg="Initializing gradient bw-prop of objective $k",
            kwargs...
        ) for (k, Ψ̃) in enumerate(bw_grad_states)
    ]
    tau_grads::Vector{Matrix{ComplexF64}} =
        [zeros(ComplexF64, length(controls), length(tlist) - 1) for _ in objectives]
    GrapeWrk(
        objectives,
        adjoint_objectives,
        kwargs,
        controls,
        pulsevals,
        gradient,
        searchdirection,
        fg_count,
        result,
        bw_states,
        fw_states,
        bw_grad_states,
        tau_grads,
        G,
        TDgradG,
        gradG,
        control_derivs,
        vals_dict,
        fw_storage,
        fw_prop_wrk,
        grad_prop_wrk,
        use_threads
    )
end

function _init_objgradpropwrk(
    Ψ̃,
    tlist,
    method,
    generator;
    verbose=false,
    info_msg=nothing,
    kwargs...
)
    if verbose && !isnothing(info_msg)
        @info info_msg
    end
    initpropwrk(Ψ̃, tlist, method, generator; verbose=verbose, kwargs...)
end


text:
import POMDPs.initialstate
const IVec8 = SVector{8, Int}

@with_kw struct AODiscreteVDPTagPOMDP <: POMDP{TagState, TagAction, IVec8}
    cpomdp::VDPTagPOMDP             = VDPTagPOMDP()
    angles::Array{Float64, 1}       = range(0, stop=2*pi, length=11)[1:end-1]
    binsize::Float64                = 0.5
end
AODiscreteVDPTagPOMDP(p::VDPTagPOMDP, n_angles::Int, binsize::Float64) = AODiscreteVDPTagPOMDP(p, range(0, stop=2*pi, length=n_angles+1)[1:end-1], binsize)

@with_kw struct ADiscreteVDPTagPOMDP <: POMDP{TagState, TagAction, Vec8}
    cpomdp::VDPTagPOMDP             = VDPTagPOMDP()
    angles::Array{Float64, 1}       = range(0, stop=2*pi, length=11)[1:end-1]
end
ADiscreteVDPTagPOMDP(p::VDPTagPOMDP, n_angles::Int) = ADiscreteVDPTagPOMDP(p, range(0, stop=2*pi, length=n_angles+1)[1:end-1])

const DiscreteVDPTagProblem = Union{AODiscreteVDPTagPOMDP, ADiscreteVDPTagPOMDP}

cproblem(p::DiscreteVDPTagProblem) = p.cpomdp::VDPTagPOMDP
mdp(p::DiscreteVDPTagProblem) = mdp(cproblem(p))::VDPTagMDP

convert_s(::Type{T}, x::T, p) where T = x
convert_a(::Type{T}, x::T, p) where T = x
convert_o(::Type{T}, x::T, p) where T = x

# observation
function convert_o(::Type{IVec8}, o::Vec8, p::AODiscreteVDPTagPOMDP)
    return floor.(Int, (o./p.binsize)::Vec8)::IVec8
end
# convert_o(::Type{Vec8}, o::Int, p::DiscreteVDPTagProblem) = (o-0.5)*2*pi/p.n_obs_angles

n_states(p::DiscreteVDPTagProblem) = Inf
n_actions(p::DiscreteVDPTagProblem) = 2*length(p.angles)
POMDPs.discount(p::DiscreteVDPTagProblem) = discount(cproblem(p))
POMDPs.isterminal(p::DiscreteVDPTagProblem, s::TagState) = mdp(p).tag_terminate && norm(s.agent-s.target) < mdp(p).tag_radius
POMDPs.actions(p::DiscreteVDPTagProblem) = [TagAction(look, angle) for look in [false, true] for angle in p.angles]
POMDPs.actionindex(p::DiscreteVDPTagProblem, a::TagAction) = a.look * length(p.angles) + findfirst(x->x==a.angle, p.angles)

POMDPs.transition(p::DiscreteVDPTagProblem, s::TagState, a::TagAction) = transition(cproblem(p), s, a)
POMDPs.initialstate(p::DiscreteVDPTagProblem) = VDPInitDist()
POMDPs.reward(p::DiscreteVDPTagProblem, s::TagState, a::TagAction, sp::TagState) = reward(cproblem(p), s, a, sp)

POMDPs.observation(p::ADiscreteVDPTagPOMDP, a::TagAction, sp::TagState) = observation(cproblem(p), a, sp)

struct DiscreteBeamDist
    beam_dist::BeamDist
    pomdp::POMDP
end

POMDPs.observation(p::AODiscreteVDPTagPOMDP, a::TagAction, sp::TagState) = DiscreteBeamDist(observation(cproblem(p), a, sp), p)
rand(rng::AbstractRNG, d::DiscreteBeamDist) = convert_o(IVec8, rand(rng, d.beam_dist), d.pomdp)

function POMDPs.pdf(d::DiscreteBeamDist, o::IVec8)
    p = 1.0
    lower = o .* d.pomdp.binsize
    upper = (o .+ 1) .* d.pomdp.binsize
    d = d.beam_dist
    for i in 1:length(o)
        if i == d.abeam
            p *= cdf(d.an, upper[i]) - cdf(d.an, lower[i])
        else
            p *= cdf(d.n, upper[i]) - cdf(d.n, lower[i])
        end
    end
    return p
end