Using quadgk inside a custom objective

[unary-code]

When I run this code:

jl.eval("""
x1 = Node(;feature=1)

tree = (((((x1 + x1) - x1) / ((x1 * x1) + (x1 - 0.101250745))) / 0.3167079) + (0.101250745 / -0.61750597))

tree = Node{Float32}(tree)

global integral = 1
err = -3
try
    my_tuple = quadgk(x -> (eval_tree_array(tree, reshape([Float32(x)], 1, 1), options))[1][1], 0, 1, rtol=0.001)
    global integral = my_tuple[1]
    global err = my_tuple[2]
    println("integral so the result of calling quadgk the first time is", integral)
catch
    println("A error when calling quadgk the first time")
end



try
    my_tuple_new = quadgk(x -> (((((x + x) - x) / ((x * x) + (x - 0.101250745))) / 0.3167079) + (0.101250745 / -0.61750597)), 0, 1, rtol=0.001)
    println("integral so the result of calling quadgk the first time is", my_tuple_new[1])

catch
    println("A error when calling quadgk the second time")
end





""")

this is the output:

integral so the result of calling quadgk the first time is0.005200976747003625
A error when calling quadgk the second time

I want a error to happen because I know that the definite integral of this math function from x=0 to x=1 is infinity and/or undefined.
Unfortunately, in my Julia objective function, I am currently passing eval_tree_array into the quadgk function, which means that when I do "model.fit(X, y)", my Julia objective function actually finds a finite value (in this case, 0.0052009), so my Julia objective function actually returns a finite value for the loss.

My guess for why this difference happens is that eval_tree_array is black-box from the perspective of the quadgk function.
In other words, the quadgk function does not understand what the symbolic form of eval_tree_array in terms of x, which means that quadgk has to do a numerical Riemann-sum approximation when you pass eval_tree_array into it, so the output of quadgk is less accurate when you pass eval_tree_array into the quadgk function, then when you call quadgk function with a direct symbolic expression in terms of "x".

Any ways to get around this, so that the Julia objective function can be accurate (can actually throw a error when it calls quadgk, like in the second time I called quadgk above)?
Perhaps I could pass the string output of "string_tree(tree, options)" into quadgk in the Julia objective function?

[MilesCranmer]

Looking now. General advice is to avoid wrapping things with a catch-all try-catch otherwise it might just skip an error you didn't expect. It's better to only catch the errors you want to, like

try
    error("...")
catch e
    isa(e, ErrorException) || rethrow(e)
end

[MilesCranmer]

Part of the issue is that you are avoiding checking the flag for eval_tree_array:

eval_tree_array(tree, reshape([Float32(x)], 1, 1), options))[1][1]

The return of eval_tree_array contains two things. The second thing it returns, which is being avoided here with the [1], is the completion flag. The completion flag is false if an infinity is encountered during evaluation.

What you actually want is something like

my_tuple = quadgk(
    x -> begin
       output, completion = (eval_tree_array(tree, reshape([Float32(x)], 1, 1), options))
       if !completion
           error("Integration failed.")
       end
       output[1]
    end,
    0,
    1,
    rtol=0.001
)

[MilesCranmer]

It actually looks like there is no way for QuadGK.jl to detect singularities: https://juliamath.github.io/QuadGK.jl/stable/quadgk-examples/. It will integrate right through them.

The way it recommends to get around is to look at the number of steps it takes. So for example, rather than setting rtol, you could set the maxevals, and just assume there's a singularity if it needs more evals than that:

function integrate_tree(tree, options)
    result, result_err, actual_evals = try
        quadgk_count(
            x -> begin
               output, completion = eval_tree_array(tree, reshape([Float32(x)], 1, 1), options)
               if !completion
                   error("Integration failed.")
               end
               output[1]
            end,
            0,
            1; 
            maxevals=1000
        )
    catch e
        isa(e, ErrorException) || rethrow(e)
        return L(Inf)
    end
    
    if actual_evals >= 1000
        return L(Inf)
    end
    return result # `result` stores the integration value
end

make sense?

For other issues about QuadGK.jl you might want to ask on their issues page as it is reaching the limits of my knowledge.

Cheers,
Miles

[MilesCranmer]

Works for me:

julia> integrate_tree(Node(; val=1.0f0), options)
1.0

julia> integrate_tree(x1, options)
0.5

julia> integrate_tree(2 * x1, options)
1.0

julia> integrate_tree(2 * x1 / (x1 * x1), options)
Inf32