diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 5c9331c9..05ea4678 100644
--- a/benchmark/benchmarks.jl
+++ b/benchmark/benchmarks.jl
@@ -74,11 +74,18 @@ function benchmark_evaluation()
                 extra_kws...
             )
             suite[T]["evaluation$(extra_key)"] = @benchmarkable(
-                [eval_tree_array(tree, X, $operators; turbo=$turbo, $extra_kws...) for tree in trees],
+                [eval_tree_array(tree, X, $operators; kws...) for tree in trees],
                 setup=(
                     X=randn(MersenneTwister(0), $T, 5, $n);
                     treesize=20;
                     ntrees=100;
+                    kws=$(
+                        if @isdefined(EvalOptions)
+                            (; eval_options=EvalOptions(; turbo=turbo, extra_kws...))
+                        else
+                            (; turbo, extra_kws...)
+                        end
+                    );
                     trees=[gen_random_tree_fixed_size(treesize, $operators, 5, $T) for _ in 1:ntrees]
                 )
             )
diff --git a/docs/src/eval.md b/docs/src/eval.md
index 82adec7b..370cda80 100644
--- a/docs/src/eval.md
+++ b/docs/src/eval.md
@@ -6,14 +6,18 @@ Given an expression tree specified with a `Node` type, you may evaluate the expr
 over an array of data with the following command:
 
 ```@docs
-eval_tree_array(tree::Node{T}, cX::AbstractMatrix{T}, operators::OperatorEnum) where {T<:Number}
+eval_tree_array(
+    tree::AbstractExpressionNode{T},
+    cX::AbstractMatrix{T},
+    operators::OperatorEnum;
+    eval_options::Union{EvalOptions,Nothing}=nothing,
+) where {T}
 ```
 
-Assuming you are only using a single `OperatorEnum`, you can also use
-the following shorthand by using the expression as a function:
+You can also use the following shorthand by using the expression as a function:
 
 ```
-    (tree::AbstractExpressionNode)(X::AbstractMatrix{T}, operators::OperatorEnum; turbo::Union{Bool,Val}=false, bumper::Union{Bool,Val}=Val(false))
+    (tree::AbstractExpressionNode)(X, operators::OperatorEnum; kws...)
 
 Evaluate a binary tree (equation) over a given input data matrix. The
 operators contain all of the operators used. This function fuses doublets
@@ -23,8 +27,7 @@ and triplets of operations for lower memory usage.
 - `tree::AbstractExpressionNode`: The root node of the tree to evaluate.
 - `cX::AbstractMatrix{T}`: The input data to evaluate the tree on.
 - `operators::OperatorEnum`: The operators used in the tree.
-- `turbo::Union{Bool,Val}`: Use LoopVectorization.jl for faster evaluation.
-- `bumper::Union{Bool,Val}`: Use Bumper.jl for faster evaluation.
+- `kws...`: Passed to [`eval_tree_array`](@ref).
 
 # Returns
 - `output::AbstractVector{T}`: the result, which is a 1D array.
@@ -53,6 +56,14 @@ It also re-defines `print`, `show`, and the various operators, to work with the
     Thus, if you define an expression with one `OperatorEnum`, and then try to
     evaluate it or print it with a different `OperatorEnum`, you will get undefined behavior!
 
+    For safer behavior, you should use [`Expression`](@ref) objects.
+
+Evaluation options are specified using `EvalOptions`:
+
+```@docs
+EvalOptions
+```
+
 You can also work with arbitrary types, by defining a `GenericOperatorEnum` instead.
 The notation is the same for `eval_tree_array`, though it will return `nothing`
 when it can't find a method, and not do any NaN checks:
diff --git a/ext/DynamicExpressionsBumperExt.jl b/ext/DynamicExpressionsBumperExt.jl
index 0b0c99f3..6e99927b 100644
--- a/ext/DynamicExpressionsBumperExt.jl
+++ b/ext/DynamicExpressionsBumperExt.jl
@@ -2,7 +2,7 @@ module DynamicExpressionsBumperExt
 
 using Bumper: @no_escape, @alloc
 using DynamicExpressions:
-    OperatorEnum, AbstractExpressionNode, tree_mapreduce, is_valid_array
+    OperatorEnum, AbstractExpressionNode, tree_mapreduce, is_valid_array, EvalOptions
 using DynamicExpressions.UtilsModule: ResultOk, counttuple
 
 import DynamicExpressions.ExtensionInterfaceModule:
@@ -12,8 +12,8 @@ function bumper_eval_tree_array(
     tree::AbstractExpressionNode{T},
     cX::AbstractMatrix{T},
     operators::OperatorEnum,
-    ::Val{turbo},
-) where {T,turbo}
+    eval_options::EvalOptions{turbo,true,early_exit},
+) where {T,turbo,early_exit}
     result = similar(cX, axes(cX, 2))
     n = size(cX, 2)
     all_ok = Ref(false)
@@ -26,7 +26,7 @@ function bumper_eval_tree_array(
                 ok = if leaf_node.constant
                     v = leaf_node.val
                     ar .= v
-                    isfinite(v)
+                    early_exit ? isfinite(v) : true
                 else
                     ar .= view(cX, leaf_node.feature, :)
                     true
@@ -38,7 +38,7 @@ function bumper_eval_tree_array(
             # In the evaluation kernel, we combine the branch nodes
             # with the arrays created by the leaf nodes:
             ((args::Vararg{Any,M}) where {M}) ->
-                dispatch_kerns!(operators, args..., Val(turbo)),
+                dispatch_kerns!(operators, args..., eval_options),
             tree;
             break_sharing=Val(true),
         )
@@ -49,55 +49,61 @@ function bumper_eval_tree_array(
     return (result, all_ok[])
 end
 
-function dispatch_kerns!(operators, branch_node, cumulator, ::Val{turbo}) where {turbo}
+function dispatch_kerns!(
+    operators, branch_node, cumulator, eval_options::EvalOptions{<:Any,true,early_exit}
+) where {early_exit}
     cumulator.ok || return cumulator
 
-    out = dispatch_kern1!(operators.unaops, branch_node.op, cumulator.x, Val(turbo))
-    return ResultOk(out, is_valid_array(out))
+    out = dispatch_kern1!(operators.unaops, branch_node.op, cumulator.x, eval_options)
+    return ResultOk(out, early_exit ? is_valid_array(out) : true)
 end
 function dispatch_kerns!(
-    operators, branch_node, cumulator1, cumulator2, ::Val{turbo}
-) where {turbo}
+    operators,
+    branch_node,
+    cumulator1,
+    cumulator2,
+    eval_options::EvalOptions{<:Any,true,early_exit},
+) where {early_exit}
     cumulator1.ok || return cumulator1
     cumulator2.ok || return cumulator2
 
     out = dispatch_kern2!(
-        operators.binops, branch_node.op, cumulator1.x, cumulator2.x, Val(turbo)
+        operators.binops, branch_node.op, cumulator1.x, cumulator2.x, eval_options
     )
-    return ResultOk(out, is_valid_array(out))
+    return ResultOk(out, early_exit ? is_valid_array(out) : true)
 end
 
-@generated function dispatch_kern1!(unaops, op_idx, cumulator, ::Val{turbo}) where {turbo}
+@generated function dispatch_kern1!(unaops, op_idx, cumulator, eval_options::EvalOptions)
     nuna = counttuple(unaops)
     quote
         Base.@nif(
             $nuna,
             i -> i == op_idx,
             i -> let op = unaops[i]
-                return bumper_kern1!(op, cumulator, Val(turbo))
+                return bumper_kern1!(op, cumulator, eval_options)
             end,
         )
     end
 end
 @generated function dispatch_kern2!(
-    binops, op_idx, cumulator1, cumulator2, ::Val{turbo}
-) where {turbo}
+    binops, op_idx, cumulator1, cumulator2, eval_options::EvalOptions
+)
     nbin = counttuple(binops)
     quote
         Base.@nif(
             $nbin,
             i -> i == op_idx,
             i -> let op = binops[i]
-                return bumper_kern2!(op, cumulator1, cumulator2, Val(turbo))
+                return bumper_kern2!(op, cumulator1, cumulator2, eval_options)
             end,
         )
     end
 end
-function bumper_kern1!(op::F, cumulator, ::Val{false}) where {F}
+function bumper_kern1!(op::F, cumulator, ::EvalOptions{false,true}) where {F}
     @. cumulator = op(cumulator)
     return cumulator
 end
-function bumper_kern2!(op::F, cumulator1, cumulator2, ::Val{false}) where {F}
+function bumper_kern2!(op::F, cumulator1, cumulator2, ::EvalOptions{false,true}) where {F}
     @. cumulator1 = op(cumulator1, cumulator2)
     return cumulator1
 end
diff --git a/ext/DynamicExpressionsLoopVectorizationExt.jl b/ext/DynamicExpressionsLoopVectorizationExt.jl
index ec158320..7edbd704 100644
--- a/ext/DynamicExpressionsLoopVectorizationExt.jl
+++ b/ext/DynamicExpressionsLoopVectorizationExt.jl
@@ -3,7 +3,7 @@ module DynamicExpressionsLoopVectorizationExt
 using LoopVectorization: @turbo
 using DynamicExpressions: AbstractExpressionNode
 using DynamicExpressions.UtilsModule: ResultOk, fill_similar
-using DynamicExpressions.EvaluateModule: @return_on_check
+using DynamicExpressions.EvaluateModule: @return_on_nonfinite_val, EvalOptions
 import DynamicExpressions.EvaluateModule:
     deg1_eval,
     deg2_eval,
@@ -18,7 +18,10 @@ import DynamicExpressions.ExtensionInterfaceModule:
 _is_loopvectorization_loaded(::Int) = true
 
 function deg2_eval(
-    cumulator_l::AbstractVector{T}, cumulator_r::AbstractVector{T}, op::F, ::Val{true}
+    cumulator_l::AbstractVector{T},
+    cumulator_r::AbstractVector{T},
+    op::F,
+    ::EvalOptions{true},
 )::ResultOk where {T<:Number,F}
     @turbo for j in eachindex(cumulator_l)
         x = op(cumulator_l[j], cumulator_r[j])
@@ -28,7 +31,7 @@ function deg2_eval(
 end
 
 function deg1_eval(
-    cumulator::AbstractVector{T}, op::F, ::Val{true}
+    cumulator::AbstractVector{T}, op::F, ::EvalOptions{true}
 )::ResultOk where {T<:Number,F}
     @turbo for j in eachindex(cumulator)
         x = op(cumulator[j])
@@ -38,21 +41,25 @@ function deg1_eval(
 end
 
 function deg1_l2_ll0_lr0_eval(
-    tree::AbstractExpressionNode{T}, cX::AbstractMatrix{T}, op::F, op_l::F2, ::Val{true}
+    tree::AbstractExpressionNode{T},
+    cX::AbstractMatrix{T},
+    op::F,
+    op_l::F2,
+    eval_options::EvalOptions{true},
 ) where {T<:Number,F,F2}
     if tree.l.l.constant && tree.l.r.constant
         val_ll = tree.l.l.val
         val_lr = tree.l.r.val
-        @return_on_check val_ll cX
-        @return_on_check val_lr cX
+        @return_on_nonfinite_val(eval_options, val_ll, cX)
+        @return_on_nonfinite_val(eval_options, val_lr, cX)
         x_l = op_l(val_ll, val_lr)::T
-        @return_on_check x_l cX
+        @return_on_nonfinite_val(eval_options, x_l, cX)
         x = op(x_l)::T
-        @return_on_check x cX
+        @return_on_nonfinite_val(eval_options, x, cX)
         return ResultOk(fill_similar(x, cX, axes(cX, 2)), true)
     elseif tree.l.l.constant
         val_ll = tree.l.l.val
-        @return_on_check val_ll cX
+        @return_on_nonfinite_val(eval_options, val_ll, cX)
         feature_lr = tree.l.r.feature
         cumulator = similar(cX, axes(cX, 2))
         @turbo for j in axes(cX, 2)
@@ -64,7 +71,7 @@ function deg1_l2_ll0_lr0_eval(
     elseif tree.l.r.constant
         feature_ll = tree.l.l.feature
         val_lr = tree.l.r.val
-        @return_on_check val_lr cX
+        @return_on_nonfinite_val(eval_options, val_lr, cX)
         cumulator = similar(cX, axes(cX, 2))
         @turbo for j in axes(cX, 2)
             x_l = op_l(cX[feature_ll, j], val_lr)
@@ -86,15 +93,19 @@ function deg1_l2_ll0_lr0_eval(
 end
 
 function deg1_l1_ll0_eval(
-    tree::AbstractExpressionNode{T}, cX::AbstractMatrix{T}, op::F, op_l::F2, ::Val{true}
+    tree::AbstractExpressionNode{T},
+    cX::AbstractMatrix{T},
+    op::F,
+    op_l::F2,
+    eval_options::EvalOptions{true},
 ) where {T<:Number,F,F2}
     if tree.l.l.constant
         val_ll = tree.l.l.val
-        @return_on_check val_ll cX
+        @return_on_nonfinite_val(eval_options, val_ll, cX)
         x_l = op_l(val_ll)::T
-        @return_on_check x_l cX
+        @return_on_nonfinite_val(eval_options, x_l, cX)
         x = op(x_l)::T
-        @return_on_check x cX
+        @return_on_nonfinite_val(eval_options, x, cX)
         return ResultOk(fill_similar(x, cX, axes(cX, 2)), true)
     else
         feature_ll = tree.l.l.feature
@@ -109,20 +120,23 @@ function deg1_l1_ll0_eval(
 end
 
 function deg2_l0_r0_eval(
-    tree::AbstractExpressionNode{T}, cX::AbstractMatrix{T}, op::F, ::Val{true}
+    tree::AbstractExpressionNode{T},
+    cX::AbstractMatrix{T},
+    op::F,
+    eval_options::EvalOptions{true},
 ) where {T<:Number,F}
     if tree.l.constant && tree.r.constant
         val_l = tree.l.val
-        @return_on_check val_l cX
+        @return_on_nonfinite_val(eval_options, val_l, cX)
         val_r = tree.r.val
-        @return_on_check val_r cX
+        @return_on_nonfinite_val(eval_options, val_r, cX)
         x = op(val_l, val_r)::T
-        @return_on_check x cX
+        @return_on_nonfinite_val(eval_options, x, cX)
         return ResultOk(fill_similar(x, cX, axes(cX, 2)), true)
     elseif tree.l.constant
         cumulator = similar(cX, axes(cX, 2))
         val_l = tree.l.val
-        @return_on_check val_l cX
+        @return_on_nonfinite_val(eval_options, val_l, cX)
         feature_r = tree.r.feature
         @turbo for j in axes(cX, 2)
             x = op(val_l, cX[feature_r, j])
@@ -133,7 +147,7 @@ function deg2_l0_r0_eval(
         cumulator = similar(cX, axes(cX, 2))
         feature_l = tree.l.feature
         val_r = tree.r.val
-        @return_on_check val_r cX
+        @return_on_nonfinite_val(eval_options, val_r, cX)
         @turbo for j in axes(cX, 2)
             x = op(cX[feature_l, j], val_r)
             cumulator[j] = x
@@ -157,11 +171,11 @@ function deg2_l0_eval(
     cumulator::AbstractVector{T},
     cX::AbstractArray{T},
     op::F,
-    ::Val{true},
+    eval_options::EvalOptions{true},
 ) where {T<:Number,F}
     if tree.l.constant
         val = tree.l.val
-        @return_on_check val cX
+        @return_on_nonfinite_val(eval_options, val, cX)
         @turbo for j in eachindex(cumulator)
             x = op(val, cumulator[j])
             cumulator[j] = x
@@ -182,11 +196,11 @@ function deg2_r0_eval(
     cumulator::AbstractVector{T},
     cX::AbstractArray{T},
     op::F,
-    ::Val{true},
+    eval_options::EvalOptions{true},
 ) where {T<:Number,F}
     if tree.r.constant
         val = tree.r.val
-        @return_on_check val cX
+        @return_on_nonfinite_val(eval_options, val, cX)
         @turbo for j in eachindex(cumulator)
             x = op(cumulator[j], val)
             cumulator[j] = x
@@ -203,11 +217,15 @@ function deg2_r0_eval(
 end
 
 ## Interface with Bumper.jl
-function bumper_kern1!(op::F, cumulator, ::Val{true}) where {F}
+function bumper_kern1!(
+    op::F, cumulator, ::EvalOptions{true,true,early_exit}
+) where {F,early_exit}
     @turbo @. cumulator = op(cumulator)
     return cumulator
 end
-function bumper_kern2!(op::F, cumulator1, cumulator2, ::Val{true}) where {F}
+function bumper_kern2!(
+    op::F, cumulator1, cumulator2, ::EvalOptions{true,true,early_exit}
+) where {F,early_exit}
     @turbo @. cumulator1 = op(cumulator1, cumulator2)
     return cumulator1
 end
diff --git a/src/DynamicExpressions.jl b/src/DynamicExpressions.jl
index dbe80e2b..25ff39b0 100644
--- a/src/DynamicExpressions.jl
+++ b/src/DynamicExpressions.jl
@@ -70,7 +70,8 @@ import .NodeModule:
 @reexport import .OperatorEnumModule: AbstractOperatorEnum
 @reexport import .OperatorEnumConstructionModule:
     OperatorEnum, GenericOperatorEnum, @extend_operators, set_default_variable_names!
-@reexport import .EvaluateModule: eval_tree_array, differentiable_eval_tree_array
+@reexport import .EvaluateModule:
+    eval_tree_array, differentiable_eval_tree_array, EvalOptions
 @reexport import .EvaluateDerivativeModule: eval_diff_tree_array, eval_grad_tree_array
 @reexport import .ChainRulesModule: NodeTangent, extract_gradient
 @reexport import .SimplifyModule: combine_operators, simplify_tree!
diff --git a/src/Evaluate.jl b/src/Evaluate.jl
index 88999e62..670a983b 100644
--- a/src/Evaluate.jl
+++ b/src/Evaluate.jl
@@ -1,6 +1,6 @@
 module EvaluateModule
 
-using DispatchDoctor: @unstable
+using DispatchDoctor: @stable, @unstable
 
 import ..NodeModule: AbstractExpressionNode, constructorof
 import ..StringsModule: string_tree
@@ -12,24 +12,90 @@ import ..ValueInterfaceModule: is_valid, is_valid_array
 
 const OPERATOR_LIMIT_BEFORE_SLOWDOWN = 15
 
-macro return_on_check(val, X)
+macro return_on_nonfinite_val(eval_options, val, X)
     :(
-        if !is_valid($(esc(val)))
+        if $(esc(eval_options)).early_exit isa Val{true} && !is_valid($(esc(val)))
             return $(ResultOk)(similar($(esc(X)), axes($(esc(X)), 2)), false)
         end
     )
 end
 
-macro return_on_nonfinite_array(array)
+macro return_on_nonfinite_array(eval_options, array)
     :(
-        if !is_valid_array($(esc(array)))
+        if $(esc(eval_options)).early_exit isa Val{true} && !is_valid_array($(esc(array)))
             return $(ResultOk)($(esc(array)), false)
         end
     )
 end
 
 """
-    eval_tree_array(tree::AbstractExpressionNode, cX::AbstractMatrix{T}, operators::OperatorEnum; turbo::Union{Bool,Val}=Val(false), bumper::Union{Bool,Val}=Val(false))
+    EvalOptions{T,B,E}
+
+This holds options for expression evaluation, such as evaluation backend.
+
+# Fields
+
+- `turbo::Val{T}`: If `Val{true}`, use LoopVectorization.jl for faster
+    evaluation.
+- `bumper::Val{B}`: If `Val{true}, use Bumper.jl for faster evaluation.
+- `early_exit::Val{E}`: If `Val{true}`, any element of any step becoming
+    `NaN` or `Inf` will terminate the computation and the whole buffer will be
+    returned with `NaN`s. This makes sure that expressions with singularities
+    don't wast compute cycles. Setting `Val{false}` will continue the computation
+    as usual and thus result in `NaN`s only in the elements that actually have
+    `NaN`s.
+"""
+struct EvalOptions{T,B,E}
+    turbo::Val{T}
+    bumper::Val{B}
+    early_exit::Val{E}
+end
+
+@stable(
+    default_mode = "disable",
+    default_union_limit = 2,
+    @inline _to_bool_val(x::Bool) = x ? Val(true) : Val(false)
+)
+@inline _to_bool_val(x::Val{T}) where {T} = Val(T::Bool)
+
+@unstable function EvalOptions(;
+    turbo::Union{Bool,Val}=Val(false),
+    bumper::Union{Bool,Val}=Val(false),
+    early_exit::Union{Bool,Val}=Val(true),
+)
+    return EvalOptions(_to_bool_val(turbo), _to_bool_val(bumper), _to_bool_val(early_exit))
+end
+
+@unstable function _process_deprecated_kws(eval_options, deprecated_kws)
+    turbo = get(deprecated_kws, :turbo, nothing)
+    bumper = get(deprecated_kws, :bumper, nothing)
+    if any(Base.Fix2(∉, (:turbo, :bumper)), keys(deprecated_kws))
+        throw(ArgumentError("Invalid keyword argument(s): $(keys(deprecated_kws))"))
+    end
+    if !isempty(deprecated_kws)
+        @assert eval_options === nothing "Cannot use both `eval_options` and deprecated flags `turbo` and `bumper`."
+        Base.depwarn(
+            "The `turbo` and `bumper` keyword arguments are deprecated. Please use `eval_options` instead.",
+            :eval_tree_array,
+        )
+    end
+    if eval_options !== nothing
+        return eval_options
+    else
+        return EvalOptions(;
+            turbo=turbo === nothing ? Val(false) : turbo,
+            bumper=bumper === nothing ? Val(false) : bumper,
+        )
+    end
+end
+
+"""
+    eval_tree_array(
+        tree::AbstractExpressionNode{T},
+        cX::AbstractMatrix{T},
+        operators::OperatorEnum;
+        eval_options::Union{EvalOptions,Nothing}=nothing,
+    ) where {T}
 
 Evaluate a binary tree (equation) over a given input data matrix. The
 operators contain all of the operators used. This function fuses doublets
@@ -39,8 +105,9 @@ and triplets of operations for lower memory usage.
 - `tree::AbstractExpressionNode`: The root node of the tree to evaluate.
 - `cX::AbstractMatrix{T}`: The input data to evaluate the tree on.
 - `operators::OperatorEnum`: The operators used in the tree.
-- `turbo::Union{Bool,Val}`: Use LoopVectorization.jl for faster evaluation.
-- `bumper::Union{Bool,Val}`: Use Bumper.jl for faster evaluation.
+- `eval_options::Union{EvalOptions,Nothing}`: See [`EvalOptions`](@ref) for documentation
+    on the different evaluation modes.
+
 
 # Returns
 - `(output, complete)::Tuple{AbstractVector{T}, Bool}`: the result,
@@ -68,29 +135,32 @@ function eval_tree_array(
     tree::AbstractExpressionNode{T},
     cX::AbstractMatrix{T},
     operators::OperatorEnum;
-    turbo::Union{Bool,Val}=Val(false),
-    bumper::Union{Bool,Val}=Val(false),
+    eval_options::Union{EvalOptions,Nothing}=nothing,
+    _deprecated_kws...,
 ) where {T}
-    v_turbo = isa(turbo, Val) ? turbo : (turbo ? Val(true) : Val(false))
-    v_bumper = isa(bumper, Val) ? bumper : (bumper ? Val(true) : Val(false))
-    if v_turbo isa Val{true} || v_bumper isa Val{true}
+    _eval_options = _process_deprecated_kws(eval_options, _deprecated_kws)
+    if _eval_options.turbo isa Val{true} || _eval_options.bumper isa Val{true}
         @assert T in (Float32, Float64)
     end
-    if v_turbo isa Val{true}
+    if _eval_options.turbo isa Val{true}
         _is_loopvectorization_loaded(0) ||
             error("Please load the LoopVectorization.jl package to use this feature.")
     end
-    if (v_turbo isa Val{true} || v_bumper isa Val{true}) && !(T <: Number)
+    if (_eval_options.turbo isa Val{true} || _eval_options.bumper isa Val{true}) &&
+        !(T <: Number)
         error(
             "Bumper and LoopVectorization features are only compatible with numeric element types",
         )
     end
-    if v_bumper isa Val{true}
-        return bumper_eval_tree_array(tree, cX, operators, v_turbo)
+    if _eval_options.bumper isa Val{true}
+        return bumper_eval_tree_array(tree, cX, operators, _eval_options)
     end
 
-    result = _eval_tree_array(tree, cX, operators, v_turbo)
-    return (result.x, result.ok && is_valid_array(result.x))
+    result = _eval_tree_array(tree, cX, operators, _eval_options)
+    return (
+        result.x,
+        result.ok && (_eval_options.early_exit isa Val{false} || is_valid_array(result.x)),
+    )
 end
 
 function eval_tree_array(
@@ -103,14 +173,13 @@ function eval_tree_array(
     tree::AbstractExpressionNode{T1},
     cX::AbstractMatrix{T2},
     operators::OperatorEnum;
-    turbo::Union{Bool,Val}=Val(false),
-    bumper::Union{Bool,Val}=Val(false),
+    kws...,
 ) where {T1,T2}
     T = promote_type(T1, T2)
     @warn "Warning: eval_tree_array received mixed types: tree=$(T1) and data=$(T2)."
     tree = convert(constructorof(typeof(tree)){T}, tree)
     cX = Base.Fix1(convert, T).(cX)
-    return eval_tree_array(tree, cX, operators; turbo, bumper)
+    return eval_tree_array(tree, cX, operators; kws...)
 end
 
 get_nuna(::Type{<:OperatorEnum{B,U}}) where {B,U} = counttuple(U)
@@ -120,8 +189,8 @@ function _eval_tree_array(
     tree::AbstractExpressionNode{T},
     cX::AbstractMatrix{T},
     operators::OperatorEnum,
-    ::Val{turbo},
-)::ResultOk where {T,turbo}
+    eval_options::EvalOptions,
+)::ResultOk where {T}
     # First, we see if there are only constants in the tree - meaning
     # we can just return the constant result.
     if tree.degree == 0
@@ -133,17 +202,20 @@ function _eval_tree_array(
         return ResultOk(fill_similar(const_result.x[], cX, axes(cX, 2)), true)
     elseif tree.degree == 1
         op_idx = tree.op
-        return dispatch_deg1_eval(tree, cX, op_idx, operators, Val(turbo))
+        return dispatch_deg1_eval(tree, cX, op_idx, operators, eval_options)
     else
         # TODO - add op(op2(x, y), z) and op(x, op2(y, z))
         # op(x, y), where x, y are constants or variables.
         op_idx = tree.op
-        return dispatch_deg2_eval(tree, cX, op_idx, operators, Val(turbo))
+        return dispatch_deg2_eval(tree, cX, op_idx, operators, eval_options)
     end
 end
 
 function deg2_eval(
-    cumulator_l::AbstractVector{T}, cumulator_r::AbstractVector{T}, op::F, ::Val{false}
+    cumulator_l::AbstractVector{T},
+    cumulator_r::AbstractVector{T},
+    op::F,
+    ::EvalOptions{false},
 )::ResultOk where {T,F}
     @inbounds @simd for j in eachindex(cumulator_l)
         x = op(cumulator_l[j], cumulator_r[j])::T
@@ -152,7 +224,9 @@ function deg2_eval(
     return ResultOk(cumulator_l, true)
 end
 
-function deg1_eval(cumulator::AbstractVector{T}, op::F, ::Val{false})::ResultOk where {T,F}
+function deg1_eval(
+    cumulator::AbstractVector{T}, op::F, ::EvalOptions{false}
+)::ResultOk where {T,F}
     @inbounds @simd for j in eachindex(cumulator)
         x = op(cumulator[j])::T
         cumulator[j] = x
@@ -175,20 +249,20 @@ end
     cX::AbstractMatrix{T},
     op_idx::Integer,
     operators::OperatorEnum,
-    ::Val{turbo},
-) where {T,turbo}
+    eval_options::EvalOptions,
+) where {T}
     nbin = get_nbin(operators)
     long_compilation_time = nbin > OPERATOR_LIMIT_BEFORE_SLOWDOWN
     if long_compilation_time
         return quote
-            result_l = _eval_tree_array(tree.l, cX, operators, Val(turbo))
+            result_l = _eval_tree_array(tree.l, cX, operators, eval_options)
             !result_l.ok && return result_l
-            @return_on_nonfinite_array result_l.x
-            result_r = _eval_tree_array(tree.r, cX, operators, Val(turbo))
+            @return_on_nonfinite_array(eval_options, result_l.x)
+            result_r = _eval_tree_array(tree.r, cX, operators, eval_options)
             !result_r.ok && return result_r
-            @return_on_nonfinite_array result_r.x
+            @return_on_nonfinite_array(eval_options, result_r.x)
             # op(x, y), for any x or y
-            deg2_eval(result_l.x, result_r.x, operators.binops[op_idx], Val(turbo))
+            deg2_eval(result_l.x, result_r.x, operators.binops[op_idx], eval_options)
         end
     end
     return quote
@@ -197,28 +271,28 @@ end
             i -> i == op_idx,
             i -> let op = operators.binops[i]
                 if tree.l.degree == 0 && tree.r.degree == 0
-                    deg2_l0_r0_eval(tree, cX, op, Val(turbo))
+                    deg2_l0_r0_eval(tree, cX, op, eval_options)
                 elseif tree.r.degree == 0
-                    result_l = _eval_tree_array(tree.l, cX, operators, Val(turbo))
+                    result_l = _eval_tree_array(tree.l, cX, operators, eval_options)
                     !result_l.ok && return result_l
-                    @return_on_nonfinite_array result_l.x
+                    @return_on_nonfinite_array(eval_options, result_l.x)
                     # op(x, y), where y is a constant or variable but x is not.
-                    deg2_r0_eval(tree, result_l.x, cX, op, Val(turbo))
+                    deg2_r0_eval(tree, result_l.x, cX, op, eval_options)
                 elseif tree.l.degree == 0
-                    result_r = _eval_tree_array(tree.r, cX, operators, Val(turbo))
+                    result_r = _eval_tree_array(tree.r, cX, operators, eval_options)
                     !result_r.ok && return result_r
-                    @return_on_nonfinite_array result_r.x
+                    @return_on_nonfinite_array(eval_options, result_r.x)
                     # op(x, y), where x is a constant or variable but y is not.
-                    deg2_l0_eval(tree, result_r.x, cX, op, Val(turbo))
+                    deg2_l0_eval(tree, result_r.x, cX, op, eval_options)
                 else
-                    result_l = _eval_tree_array(tree.l, cX, operators, Val(turbo))
+                    result_l = _eval_tree_array(tree.l, cX, operators, eval_options)
                     !result_l.ok && return result_l
-                    @return_on_nonfinite_array result_l.x
-                    result_r = _eval_tree_array(tree.r, cX, operators, Val(turbo))
+                    @return_on_nonfinite_array(eval_options, result_l.x)
+                    result_r = _eval_tree_array(tree.r, cX, operators, eval_options)
                     !result_r.ok && return result_r
-                    @return_on_nonfinite_array result_r.x
+                    @return_on_nonfinite_array(eval_options, result_r.x)
                     # op(x, y), for any x or y
-                    deg2_eval(result_l.x, result_r.x, op, Val(turbo))
+                    deg2_eval(result_l.x, result_r.x, op, eval_options)
                 end
             end
         )
@@ -229,16 +303,16 @@ end
     cX::AbstractMatrix{T},
     op_idx::Integer,
     operators::OperatorEnum,
-    ::Val{turbo},
-) where {T,turbo}
+    eval_options::EvalOptions,
+) where {T}
     nuna = get_nuna(operators)
     long_compilation_time = nuna > OPERATOR_LIMIT_BEFORE_SLOWDOWN
     if long_compilation_time
         return quote
-            result = _eval_tree_array(tree.l, cX, operators, Val(turbo))
+            result = _eval_tree_array(tree.l, cX, operators, eval_options)
             !result.ok && return result
-            @return_on_nonfinite_array result.x
-            deg1_eval(result.x, operators.unaops[op_idx], Val(turbo))
+            @return_on_nonfinite_array(eval_options, result.x)
+            deg1_eval(result.x, operators.unaops[op_idx], eval_options)
         end
     end
     # This @nif lets us generate an if statement over choice of operator,
@@ -252,20 +326,20 @@ end
                     # op(op2(x, y)), where x, y, z are constants or variables.
                     l_op_idx = tree.l.op
                     dispatch_deg1_l2_ll0_lr0_eval(
-                        tree, cX, op, l_op_idx, operators.binops, Val(turbo)
+                        tree, cX, op, l_op_idx, operators.binops, eval_options
                     )
                 elseif tree.l.degree == 1 && tree.l.l.degree == 0
                     # op(op2(x)), where x is a constant or variable.
                     l_op_idx = tree.l.op
                     dispatch_deg1_l1_ll0_eval(
-                        tree, cX, op, l_op_idx, operators.unaops, Val(turbo)
+                        tree, cX, op, l_op_idx, operators.unaops, eval_options
                     )
                 else
                     # op(x), for any x.
-                    result = _eval_tree_array(tree.l, cX, operators, Val(turbo))
+                    result = _eval_tree_array(tree.l, cX, operators, eval_options)
                     !result.ok && return result
-                    @return_on_nonfinite_array result.x
-                    deg1_eval(result.x, op, Val(turbo))
+                    @return_on_nonfinite_array(eval_options, result.x)
+                    deg1_eval(result.x, op, eval_options)
                 end
             end
         )
@@ -277,8 +351,8 @@ end
     op::F,
     l_op_idx::Integer,
     binops,
-    ::Val{turbo},
-) where {T,F,turbo}
+    eval_options::EvalOptions,
+) where {T,F}
     nbin = counttuple(binops)
     # (Note this is only called from dispatch_deg1_eval, which has already
     # checked for long compilation times, so we don't need to check here)
@@ -287,7 +361,7 @@ end
             $nbin,
             j -> j == l_op_idx,
             j -> let op_l = binops[j]
-                deg1_l2_ll0_lr0_eval(tree, cX, op, op_l, Val(turbo))
+                deg1_l2_ll0_lr0_eval(tree, cX, op, op_l, eval_options)
             end,
         )
     end
@@ -298,36 +372,40 @@ end
     op::F,
     l_op_idx::Integer,
     unaops,
-    ::Val{turbo},
-)::ResultOk where {T,F,turbo}
+    eval_options::EvalOptions,
+)::ResultOk where {T,F}
     nuna = counttuple(unaops)
     quote
         Base.Cartesian.@nif(
             $nuna,
             j -> j == l_op_idx,
             j -> let op_l = unaops[j]
-                deg1_l1_ll0_eval(tree, cX, op, op_l, Val(turbo))
+                deg1_l1_ll0_eval(tree, cX, op, op_l, eval_options)
             end,
         )
     end
 end
 
 function deg1_l2_ll0_lr0_eval(
-    tree::AbstractExpressionNode{T}, cX::AbstractMatrix{T}, op::F, op_l::F2, ::Val{false}
+    tree::AbstractExpressionNode{T},
+    cX::AbstractMatrix{T},
+    op::F,
+    op_l::F2,
+    eval_options::EvalOptions{false,false},
 ) where {T,F,F2}
     if tree.l.l.constant && tree.l.r.constant
         val_ll = tree.l.l.val
         val_lr = tree.l.r.val
-        @return_on_check val_ll cX
-        @return_on_check val_lr cX
+        @return_on_nonfinite_val(eval_options, val_ll, cX)
+        @return_on_nonfinite_val(eval_options, val_lr, cX)
         x_l = op_l(val_ll, val_lr)::T
-        @return_on_check x_l cX
+        @return_on_nonfinite_val(eval_options, x_l, cX)
         x = op(x_l)::T
-        @return_on_check x cX
+        @return_on_nonfinite_val(eval_options, x, cX)
         return ResultOk(fill_similar(x, cX, axes(cX, 2)), true)
     elseif tree.l.l.constant
         val_ll = tree.l.l.val
-        @return_on_check val_ll cX
+        @return_on_nonfinite_val(eval_options, val_ll, cX)
         feature_lr = tree.l.r.feature
         cumulator = similar(cX, axes(cX, 2))
         @inbounds @simd for j in axes(cX, 2)
@@ -339,7 +417,7 @@ function deg1_l2_ll0_lr0_eval(
     elseif tree.l.r.constant
         feature_ll = tree.l.l.feature
         val_lr = tree.l.r.val
-        @return_on_check val_lr cX
+        @return_on_nonfinite_val(eval_options, val_lr, cX)
         cumulator = similar(cX, axes(cX, 2))
         @inbounds @simd for j in axes(cX, 2)
             x_l = op_l(cX[feature_ll, j], val_lr)::T
@@ -362,15 +440,19 @@ end
 
 # op(op2(x)) for x variable or constant
 function deg1_l1_ll0_eval(
-    tree::AbstractExpressionNode{T}, cX::AbstractMatrix{T}, op::F, op_l::F2, ::Val{false}
+    tree::AbstractExpressionNode{T},
+    cX::AbstractMatrix{T},
+    op::F,
+    op_l::F2,
+    eval_options::EvalOptions{false,false},
 ) where {T,F,F2}
     if tree.l.l.constant
         val_ll = tree.l.l.val
-        @return_on_check val_ll cX
+        @return_on_nonfinite_val(eval_options, val_ll, cX)
         x_l = op_l(val_ll)::T
-        @return_on_check x_l cX
+        @return_on_nonfinite_val(eval_options, x_l, cX)
         x = op(x_l)::T
-        @return_on_check x cX
+        @return_on_nonfinite_val(eval_options, x, cX)
         return ResultOk(fill_similar(x, cX, axes(cX, 2)), true)
     else
         feature_ll = tree.l.l.feature
@@ -386,20 +468,23 @@ end
 
 # op(x, y) for x and y variable/constant
 function deg2_l0_r0_eval(
-    tree::AbstractExpressionNode{T}, cX::AbstractMatrix{T}, op::F, ::Val{false}
+    tree::AbstractExpressionNode{T},
+    cX::AbstractMatrix{T},
+    op::F,
+    eval_options::EvalOptions{false,false},
 ) where {T,F}
     if tree.l.constant && tree.r.constant
         val_l = tree.l.val
-        @return_on_check val_l cX
+        @return_on_nonfinite_val(eval_options, val_l, cX)
         val_r = tree.r.val
-        @return_on_check val_r cX
+        @return_on_nonfinite_val(eval_options, val_r, cX)
         x = op(val_l, val_r)::T
-        @return_on_check x cX
+        @return_on_nonfinite_val(eval_options, x, cX)
         return ResultOk(fill_similar(x, cX, axes(cX, 2)), true)
     elseif tree.l.constant
         cumulator = similar(cX, axes(cX, 2))
         val_l = tree.l.val
-        @return_on_check val_l cX
+        @return_on_nonfinite_val(eval_options, val_l, cX)
         feature_r = tree.r.feature
         @inbounds @simd for j in axes(cX, 2)
             x = op(val_l, cX[feature_r, j])::T
@@ -410,7 +495,7 @@ function deg2_l0_r0_eval(
         cumulator = similar(cX, axes(cX, 2))
         feature_l = tree.l.feature
         val_r = tree.r.val
-        @return_on_check val_r cX
+        @return_on_nonfinite_val(eval_options, val_r, cX)
         @inbounds @simd for j in axes(cX, 2)
             x = op(cX[feature_l, j], val_r)::T
             cumulator[j] = x
@@ -434,11 +519,11 @@ function deg2_l0_eval(
     cumulator::AbstractVector{T},
     cX::AbstractArray{T},
     op::F,
-    ::Val{false},
+    eval_options::EvalOptions{false,false},
 ) where {T,F}
     if tree.l.constant
         val = tree.l.val
-        @return_on_check val cX
+        @return_on_nonfinite_val(eval_options, val, cX)
         @inbounds @simd for j in eachindex(cumulator)
             x = op(val, cumulator[j])::T
             cumulator[j] = x
@@ -460,11 +545,11 @@ function deg2_r0_eval(
     cumulator::AbstractVector{T},
     cX::AbstractArray{T},
     op::F,
-    ::Val{false},
+    eval_options::EvalOptions{false,false},
 ) where {T,F}
     if tree.r.constant
         val = tree.r.val
-        @return_on_check val cX
+        @return_on_nonfinite_val(eval_options, val, cX)
         @inbounds @simd for j in eachindex(cumulator)
             x = op(cumulator[j], val)::T
             cumulator[j] = x
@@ -674,12 +759,13 @@ function eval(current_node)
     tree::AbstractExpressionNode{T1},
     cX::AbstractArray{T2,N},
     operators::GenericOperatorEnum;
-    throw_errors::Bool=true,
+    throw_errors::Union{Val,Bool}=Val(true),
 ) where {T1,T2,N}
+    v_throw_errors = _to_bool_val(throw_errors)
     try
-        return _eval_tree_array_generic(tree, cX, operators, Val(true))
+        return _eval_tree_array_generic(tree, cX, operators, v_throw_errors)
     catch e
-        if !throw_errors
+        if v_throw_errors isa Val{false}
             return nothing, false
         end
         tree_s = string_tree(tree, operators)
diff --git a/src/EvaluationHelpers.jl b/src/EvaluationHelpers.jl
index 5e9a0b12..3762bcd0 100644
--- a/src/EvaluationHelpers.jl
+++ b/src/EvaluationHelpers.jl
@@ -8,7 +8,7 @@ import ..EvaluateDerivativeModule: eval_grad_tree_array
 
 # Evaluation:
 """
-    (tree::AbstractExpressionNode)(X::AbstractMatrix{T}, operators::OperatorEnum; turbo::Union{Bool,Val}=false, bumper::Union{Bool,Val}=Val(false))
+    (tree::AbstractExpressionNode)(X, operators::OperatorEnum; kws...)
 
 Evaluate a binary tree (equation) over a given input data matrix. The
 operators contain all of the operators used. This function fuses doublets
@@ -16,10 +16,10 @@ and triplets of operations for lower memory usage.
 
 # Arguments
 - `tree::AbstractExpressionNode`: The root node of the tree to evaluate.
-- `cX::AbstractMatrix{T}`: The input data to evaluate the tree on.
+- `X::AbstractMatrix{T}`: The input data to evaluate the tree on.
 - `operators::OperatorEnum`: The operators used in the tree.
-- `turbo::Union{Bool,Val}`: Use LoopVectorization.jl for faster evaluation.
-- `bumper::Union{Bool,Val}`: Use Bumper.jl for faster evaluation.
+- `kws...`: Passed to `eval_tree_array`.
+
 
 # Returns
 - `output::AbstractVector{T}`: the result, which is a 1D array.
@@ -32,18 +32,12 @@ function (tree::AbstractExpressionNode)(X, operators::OperatorEnum; kws...)
     return out
 end
 """
-    (tree::AbstractExpressionNode)(X::AbstractMatrix, operators::GenericOperatorEnum; throw_errors::Bool=true)
+    (tree::AbstractExpressionNode)(X, operators::GenericOperatorEnum; kws...)
 
 # Arguments
 - `X::AbstractArray`: The input data to evaluate the tree on.
 - `operators::GenericOperatorEnum`: The operators used in the tree.
-- `throw_errors::Bool=true`: Whether to throw errors
-    if they occur during evaluation. Otherwise,
-    MethodErrors will be caught before they happen and
-    evaluation will return `nothing`,
-    rather than throwing an error. This is useful in cases
-    where you are unsure if a particular tree is valid or not,
-    and would prefer to work with `nothing` as an output.
+- `kws...`: Passed to `eval_tree_array`.
 
 # Returns
 - `output`: the result of the evaluation.
diff --git a/src/Node.jl b/src/Node.jl
index 095a6bc5..4355aea5 100644
--- a/src/Node.jl
+++ b/src/Node.jl
@@ -142,7 +142,7 @@ be performed with this assumption, to preserve structure of the graph.
 ```julia
 julia> operators = OperatorEnum(;
            binary_operators=[+, -, *], unary_operators=[cos, sin]
-       );
+        );
 
 julia> x = GraphNode(feature=1)
 x1
diff --git a/src/precompile.jl b/src/precompile.jl
index a2967685..d16bc6b7 100644
--- a/src/precompile.jl
+++ b/src/precompile.jl
@@ -37,25 +37,33 @@ function test_all_combinations(; binary_operators, unary_operators, turbo, types
 
         # Trivial:
         for l in (x, c)
-            @ignore_domain_error eval_tree_array(l, X, operators; turbo=use_turbo)
+            @ignore_domain_error eval_tree_array(
+                l, X, operators; eval_options=EvalOptions(; turbo=use_turbo)
+            )
         end
 
         # Binary operators
         for i in eachindex(binops), l in (x, c), r in (x, c)
             tree = Node(i, l, r)
             tree = convert(Node{T}, tree)
-            @ignore_domain_error eval_tree_array(tree, X, operators; turbo=use_turbo)
+            @ignore_domain_error eval_tree_array(
+                tree, X, operators; eval_options=EvalOptions(; turbo=use_turbo)
+            )
         end
 
         # Unary operators
         for j in eachindex(unaops), k in eachindex(unaops), l in (x, c)
             tree = Node(j, l)
             tree = convert(Node{T}, tree)
-            @ignore_domain_error eval_tree_array(tree, X, operators; turbo=use_turbo)
+            @ignore_domain_error eval_tree_array(
+                tree, X, operators; eval_options=EvalOptions(; turbo=use_turbo)
+            )
 
             tree = Node(j, Node(k, l))
             tree = convert(Node{T}, tree)
-            @ignore_domain_error eval_tree_array(tree, X, operators; turbo=use_turbo)
+            @ignore_domain_error eval_tree_array(
+                tree, X, operators; eval_options=EvalOptions(; turbo=use_turbo)
+            )
         end
 
         # Both operators
@@ -67,11 +75,15 @@ function test_all_combinations(; binary_operators, unary_operators, turbo, types
 
             tree = Node(i, Node(j1, l), Node(j2, r))
             tree = convert(Node{T}, tree)
-            @ignore_domain_error eval_tree_array(tree, X, operators; turbo=use_turbo)
+            @ignore_domain_error eval_tree_array(
+                tree, X, operators; eval_options=EvalOptions(; turbo=use_turbo)
+            )
 
             tree = Node(j1, Node(i, l, r))
             tree = convert(Node{T}, tree)
-            @ignore_domain_error eval_tree_array(tree, X, operators; turbo=use_turbo)
+            @ignore_domain_error eval_tree_array(
+                tree, X, operators; eval_options=EvalOptions(; turbo=use_turbo)
+            )
         end
     end
     return nothing
diff --git a/test/test_deprecations.jl b/test/test_deprecations.jl
index 0b1af26b..45decb0f 100644
--- a/test/test_deprecations.jl
+++ b/test/test_deprecations.jl
@@ -1,6 +1,4 @@
-using DynamicExpressions
-using Test
-using Zygote
+using Test, DynamicExpressions, Zygote, LoopVectorization
 using Suppressor: @capture_err
 using DispatchDoctor: allow_unstable
 
@@ -47,6 +45,14 @@ if VERSION >= v"1.9"
     )
 end
 
+# Old usage of evaluation options
+if VERSION >= v"1.9-"
+    ex = Expression(Node{Float64}(; feature=1))
+    @test_logs (:warn, r"The `turbo` and `bumper` keyword arguments are deprecated.*") (ex(
+        randn(Float64, 1, 10), OperatorEnum(); turbo=true
+    ))
+end
+
 # Test deprecated modules
 logs = @capture_err begin
     @eval using DynamicExpressions.EquationModule
diff --git a/test/test_enzyme.jl b/test/test_enzyme.jl
index 3d33cc7c..42c08a43 100644
--- a/test/test_enzyme.jl
+++ b/test/test_enzyme.jl
@@ -54,7 +54,11 @@ X = [1.0; 1.0;;]
 d_tree = begin
     storage_tree = copy(tree)
     # Set all constants to zero:
-    Enzyme.make_zero!(storage_tree)
+    foreach(storage_tree) do node
+        if node.degree == 0 && node.constant
+            node.val = 0.0
+        end
+    end
     fetch(
         schedule(
             Task(64 * 1024^2) do
diff --git a/test/test_evaluation.jl b/test/test_evaluation.jl
index 0757488e..587fba2e 100644
--- a/test/test_evaluation.jl
+++ b/test/test_evaluation.jl
@@ -1,8 +1,7 @@
-using DynamicExpressions
-using Bumper
-using LoopVectorization
-using Random
-using Test
+#! format: off
+@testitem "Test validity of expression evaluation" begin
+using DynamicExpressions, Bumper, LoopVectorization, Random
+
 include("test_params.jl")
 include("tree_gen_utils.jl")
 
@@ -81,10 +80,19 @@ for turbo in [Val(false), Val(true)],
         end
     end
 end
+end
+#! format: on
+
+@testitem "Test specific branches of evaluation" begin
+    using DynamicExpressions, DynamicExpressions, Bumper, LoopVectorization
+    using DynamicExpressions.EvaluateModule: EvalOptions
+
+    include("test_params.jl")
 
-@testset "Test specific branches of evaluation" begin
-    for turbo in [false, true], T in [Float16, Float32, Float64, ComplexF32, ComplexF64]
-        turbo && !(T in (Float32, Float64)) && continue
+    for turbo in [Val(false), Val(true)],
+        T in [Float16, Float32, Float64, ComplexF32, ComplexF64]
+
+        turbo isa Val{true} && !(T in (Float32, Float64)) && continue
         # Test specific branches of evaluation code:
         # op(op(<constant>))
         local tree, operators
@@ -95,7 +103,7 @@ end
         @test repr(tree) == "cos(cos(3.0))"
         tree = convert(Node{T}, tree)
         truth = cos(cos(T(3.0f0)))
-        @test DynamicExpressions.EvaluateModule.deg1_l1_ll0_eval(tree, [zero(T)]', cos, cos, Val(turbo)).x[1] ≈
+        @test DynamicExpressions.EvaluateModule.deg1_l1_ll0_eval(tree, [zero(T)]', cos, cos, EvalOptions(; turbo)).x[1] ≈
             truth
 
         # op(<constant>, <constant>)
@@ -103,7 +111,7 @@ end
         @test repr(tree) == "3.0 + 4.0"
         tree = convert(Node{T}, tree)
         truth = T(3.0f0) + T(4.0f0)
-        @test DynamicExpressions.EvaluateModule.deg2_l0_r0_eval(tree, [zero(T)]', (+), Val(turbo)).x[1] ≈
+        @test DynamicExpressions.EvaluateModule.deg2_l0_r0_eval(tree, [zero(T)]', (+), EvalOptions(; turbo)).x[1] ≈
             truth
 
         # op(op(<constant>, <constant>))
@@ -111,7 +119,7 @@ end
         @test repr(tree) == "cos(3.0 + 4.0)"
         tree = convert(Node{T}, tree)
         truth = cos(T(3.0f0) + T(4.0f0))
-        @test DynamicExpressions.EvaluateModule.deg1_l2_ll0_lr0_eval(tree, [zero(T)]', cos, (+), Val(turbo)).x[1] ≈
+        @test DynamicExpressions.EvaluateModule.deg1_l2_ll0_lr0_eval(tree, [zero(T)]', cos, (+), EvalOptions(; turbo)).x[1] ≈
             truth
 
         # Test for presence of NaNs:
@@ -126,17 +134,20 @@ end
 end
 
 # Check if julia version >= 1.7:
-if VERSION >= v"1.7"
-    @testset "Test error catching for GenericOperatorEnum" begin
+@testitem "Test error catching for GenericOperatorEnum" begin
+    using DynamicExpressions
+
+    @static if VERSION >= v"1.7"
         # And, with generic operator enum, this should be an actual error:
+        @eval my_fnc(x::Real) = x
         operators = GenericOperatorEnum(;
-            binary_operators=[+, -, *, /], unary_operators=[cos, sin]
+            binary_operators=[+, -, *, /], unary_operators=[cos, sin, my_fnc]
         )
+        @extend_operators operators
         x1 = Node(Float64; feature=1)
         tree = sin(x1 / 0.0)
         X = randn(Float32, 10)
         let
-            local stack
             try
                 tree(X, operators)[1]
                 @test false
@@ -145,30 +156,34 @@ if VERSION >= v"1.7"
                 # Check that "Failed to evaluate" is in the message:
                 @test occursin("Failed to evaluate", e.msg)
                 stack = current_exceptions()
+                @test length(stack) == 2
+                @test stack[1].exception isa DomainError
             end
-            @test length(stack) == 2
-            @test stack[1].exception isa DomainError
 
             # If a method is not defined, we should get a nothing:
-            X = randn(Float32, 1, 10)
-            @test tree(X, operators; throw_errors=false) === nothing
+            X2 = randn(ComplexF64, 1, 10)
+            tree2 = my_fnc(x1)
+            @test tree2(X2, operators; throw_errors=false) === nothing
             # or a MethodError:
             try
-                tree(X, operators; throw_errors=true)
+                tree2(X2, operators; throw_errors=true)
                 @test false
             catch e
                 @test e isa ErrorException
                 @test occursin("Failed to evaluate", e.msg)
-                stack = current_exceptions()
+                stack2 = current_exceptions()
+                @test length(stack2) == 2
+                @test stack2[1].exception isa MethodError
             end
-            @test length(stack) == 2
-            # Dividing by 0 should not be an MethodError
-            # @test stack[1].exception isa MethodError
         end
     end
 end
 
-@testset "Test many operators" begin
+@testitem "Test many operators" begin
+    using DynamicExpressions
+
+    include("tree_gen_utils.jl")
+
     # Since we use `@nif` in evaluating expressions,
     # we can see if there are any issues with LARGE numbers of operators.
     num_ops = 100
@@ -211,11 +226,64 @@ end
     num_tests = 100
     n_features = 3
     for _ in 1:num_tests
-        tree = gen_random_tree_fixed_size(20, only_basic_ops_operator, n_features, Float64)
-        X = randn(Float64, n_features, 10)
-        basic_eval = tree(X, only_basic_ops_operator)
-        many_ops_eval = tree(X, many_ops_operators)
-        @test (all(isnan, basic_eval) && all(isnan, many_ops_eval)) ||
-            basic_eval ≈ many_ops_eval
+        let tree = gen_random_tree_fixed_size(
+                20, only_basic_ops_operator, n_features, Float64
+            ),
+            X = randn(Float64, n_features, 10),
+            basic_eval = tree(X, only_basic_ops_operator),
+            many_ops_eval = tree(X, many_ops_operators)
+
+            @test (all(isnan, basic_eval) && all(isnan, many_ops_eval)) ||
+                basic_eval ≈ many_ops_eval
+        end
     end
 end
+
+@testitem "Disable early exit" begin
+    using DynamicExpressions
+    using Bumper, LoopVectorization
+
+    let
+        T = Float16
+        ex = @parse_expression(
+            2 * x, binary_operators = [*], variable_names = ["x"], node_type = Node{T}
+        )
+        X = T[1.0 floatmax(T)]
+        @test all(isnan.(ex(X)))
+        @test ex(X; eval_options=EvalOptions(; early_exit=Val(false))) ≈ [2.0, Inf]
+    end
+
+    for turbo in [Val(false), Val(true)],
+        T in [Float32, Float64],
+        bumper in [Val(false), Val(true)]
+
+        ex = @parse_expression(
+            (-b - sqrt(b^2 - (4 * a) * c)) / (2 * c),
+            binary_operators = [-, *, /, ^],
+            unary_operators = [-, sqrt],
+            variable_names = ["a", "b", "c"],
+            node_type = Node{T}
+        )
+        X = T[
+            -1 -1
+            1 floatmax(T)
+            1 1
+        ]
+        @test all(isnan.(ex(X; eval_options=EvalOptions(; bumper, turbo))))
+        y = ex(X; eval_options=EvalOptions(; bumper, turbo, early_exit=Val(false)))
+        @test y[1] ≈ T(-1.618033988749895)
+        @test !isfinite(y[2])
+    end
+end
+
+@testitem "Test EvalOptions constructor" begin
+    using DynamicExpressions, LoopVectorization
+
+    @test EvalOptions(; turbo=true) isa EvalOptions{true}
+    @test EvalOptions(; turbo=Val(true)) isa EvalOptions{true}
+    @test EvalOptions(; turbo=false) isa EvalOptions{false}
+    @test EvalOptions(; turbo=Val(false)) isa EvalOptions{false}
+
+    ex = Expression(Node{Float64}(; feature=1))
+    @test_throws ArgumentError ex(randn(1, 5), OperatorEnum(); bad_arg=1)
+end
diff --git a/test/test_initial_errors.jl b/test/test_initial_errors.jl
index f5b19710..56e69cb2 100644
--- a/test/test_initial_errors.jl
+++ b/test/test_initial_errors.jl
@@ -1,4 +1,5 @@
 using DynamicExpressions
+using DynamicExpressions: EvalOptions
 using DispatchDoctor: allow_unstable
 using Test
 
@@ -39,11 +40,15 @@ if VERSION >= v"1.9"
 
     @test_throws(
         "Please load the Bumper.jl package",
-        allow_unstable(() -> tree(ones(2, 10), operators; bumper=Val(true)))
+        allow_unstable(
+            () -> tree(ones(2, 10), operators; eval_options=EvalOptions(; bumper=Val(true)))
+        )
     )
 
     @test_throws(
         "Please load the LoopVectorization.jl package",
-        allow_unstable(() -> tree(ones(2, 10), operators; turbo=Val(true)))
+        allow_unstable(
+            () -> tree(ones(2, 10), operators; eval_options=EvalOptions(; turbo=Val(true)))
+        )
     )
 end
diff --git a/test/unittest.jl b/test/unittest.jl
index 2263621f..45cb6e2e 100644
--- a/test/unittest.jl
+++ b/test/unittest.jl
@@ -46,9 +46,7 @@ end
     include("test_print.jl")
 end
 
-@testitem "Test validity of expression evaluation" begin
-    include("test_evaluation.jl")
-end
+include("test_evaluation.jl")
 
 @testitem "Test validity of integer expression evaluation" begin
     include("test_integer_evaluation.jl")