struct SimpleTauLeaping <: DiffEqBase.DEAlgorithm end

function DiffEqBase.solve(jump_prob::JumpProblem,alg::SimpleTauLeaping;
                            seed = nothing,
                            dt = error("dt is required for SimpleTauLeaping."))

    # boilerplate from SimpleTauLeaping method
    @assert isempty(jump_prob.jump_callback.continuous_callbacks) # still needs to be a regular jump
    @assert isempty(jump_prob.jump_callback.discrete_callbacks)
    prob = jump_prob.prob
    seed === nothing ? rng = Xorshifts.Xoroshiro128Plus() : rng = Xorshifts.Xoroshiro128Plus(seed)

    rj = jump_prob.regular_jump
    rate = rj.rate # rate function rate(out,u,p,t)
    numjumps = rj.numjumps # used for size information (# of jump processes)
    c = rj.c # matrix-free operator c(u_buffer, uprev, tprev, counts, p, mark)
  
    if !isnothing(rj.mark_dist) == nothing # https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/250
      error("Mark distributions are currently not supported in SimpleTauLeaping")
    end
  
    u0 = copy(prob.u0)
    du = similar(u0)
    rate_cache = zeros(float(eltype(u0)), numjumps)

    tspan = prob.tspan
    p = prob.p
  
    n = Int((tspan[2] - tspan[1])/dt) + 1
    u = Vector{typeof(prob.u0)}(undef,n)
    u[1] = u0
    t = tspan[1]:dt:tspan[2]

    # iteration variables
    counts = zero(rate_cache) # counts for each variable

    for i in 2:n # iterate over dt-slices
        uprev = u[i-1]
        tprev = t[i-1]
        rate(rate_cache,uprev,p,tprev)
        rate_cache .*= dt # multiply by the width of the time interval
        counts .= pois_rand.((rng,), rate_cache) # set counts to the poisson arrivals with our given rates
        c(du, uprev, p, tprev, counts, mark)
        u[i] = du + uprev
    end

    sol = DiffEqBase.build_solution(prob,alg,t,u,
                                    calculate_error = false,
                                    interp = DiffEqBase.ConstantInterpolation(t,u))
end

export SimpleTauLeaping
