Custom objective questions

[unary-code]

EDIT: Resolved. I figured out the reason for the difference.
In my Julia helper's code, length(prediction) was 100 (correct, because I am trying to divide by the number of points, was 100).
While in my model's objective function, length(prediction) was 1.
So, in the model's objective function, I just changed length(prediction) to sizeof(prediction)[1].

Hi,

When I run the following objective function:

objective = """
import Pkg
Pkg.add("Roots")
Pkg.add("ForwardDiff")
Pkg.add("QuadGK")

using Roots
using ForwardDiff
using QuadGK
function eval_loss(tree, dataset::Dataset{T,L}, options)::L where {T,L}
    MAX_THRESHOLD_ALLOWED = 10000
    

    prediction, flag = eval_tree_array(tree, dataset.X, options)
    !flag && return L(Inf)
    
    num_points = length(prediction)
    
    for i in 1:num_points
        cur_expr_val = prediction[i]
        if (cur_expr_val > MAX_THRESHOLD_ALLOWED)
            return L(Inf)
        end
        if (cur_expr_val < 0)
            return L(Inf)
        end
    end
    
    
    
    
    # BELOW SECTION IS TO CHECK FOR IF THE GUESSED FUNCTION TAKES ON ANY POSITIVE INFINITY VALUES
    # OR NEGATIVE INFINITY VALUES ON THE DOMAIN X=[0, 1].
    
    
    # NOTE: reshape([x], 1, 1) should give same dimensions
    # as reshape([x;], 1, 1).
    # f(x) = x -> 1/(eval_tree_array(tree, reshape([Float32(x)], 1, 1), options))
    f(x) = 1 / ((eval_tree_array(tree, reshape([Float32(x)], 1, 1), options))[1][1])
    
    # fp(x) = ForwardDiff.derivative(f,float(x))
    fp(x) = ( -1 / ((eval_tree_array(tree, reshape([Float32(x)], 1, 1), options))[1][1])^2 * (eval_diff_tree_array(tree, reshape([Float32(x)], 1, 1), options, 1))[2][1])
    
    global found_a_recip_root = true
    global recip_roots = -3
    try
        global recip_roots = find_zero((f, fp), 0.3, Roots.Newton(), maxevals=100)
        #println("recip_roots is", recip_roots)
    catch
        global found_a_recip_root = false
        # "Convergence Failed"
    end
    
    # If the function f(x) represented by the tree "tree"
    # has a value f(x) close to positive infinity or negative infinity
    # when x is in [0, 1], then return L(Inf).
    if found_a_recip_root
        if (recip_roots >= 0 && recip_roots <= 1)
            return L(Inf)
        end
    end
    
    
    # BELOW SECTION IS TO CHECK FOR IF THE GUESSED FUNCTION TAKES ON ANY NEGATIVE VALUES OR ZERO VALUES
    # ON THE DOMAIN X=[0, 1].
    f(x) = ((eval_tree_array(tree, reshape([Float32(x)], 1, 1), options))[1][1])
    
    fp(x) = ((eval_diff_tree_array(tree, reshape([Float32(x)], 1, 1), options, 1))[2][1])
    
    global found_a_root = true
    global roots = -3
    try
        global roots = find_zero((f, fp), 0.3, Roots.Newton(), maxevals=100)
        #println("roots is", roots)
    catch
        global found_a_root = false
        # "Convergence Failed"
    end
    
    if found_a_root
        if (roots >= 0 && roots <= 1)
            println("found a root so returned L(Inf) and string rep of the tree was")
            println(string_tree(tree, options))
            return L(Inf)
        end
    end
    
    
    
    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.01)
        global integral = my_tuple[1]
        global err = my_tuple[2]
    catch
        return L(Inf)
    end
    
    
    
    norm_constant = integral
    if (integral <= 0)
        println("integral was <=0")
        return L(Inf)
    end
    
    if (isinf(integral))
        #println("integral was Inf")
        return L(Inf)
    end
    
    actual_probs = (prediction) / norm_constant
    
    # length(actual_probs) equals length(prediction)
    
    # We really care about the EXPECTED VALUE OF L, where L = -1 * sum(log.(actual_probs)) / (length(prediction))
    # given that dataset.X is generated by the true probability distribution.
    # E(L) is really equal to the definite integral from x=0 to x=1 of ln(guessed_f(x)) * f_true(x) dx.
    # As num_points goes to +infinity, the variance of L becomes 0.
    
    final_loss = exp(-1 * sum(log.(actual_probs)) / (length(prediction)))
    #println("final_loss is ", final_loss)
    return final_loss
end
"""

and then I define the PySR model as:

model = PySRRegressor(
    niterations=20,  # < Increase me for better results
    binary_operators=["+", "*", "-", "/"],
    full_objective=objective,
    enable_autodiff=True,
    early_stop_condition = "f(loss, complexity) = (loss < 0.1)"
)

I define the dataset in Python as:

import numpy as np
X = np.random.rand(100, 1)
print(X.shape)
num_points = X.shape[0] # num_points is 100
y = 1000 * np.ones(num_points)

Then I run

model.fit(X, y)

(I do all this in a Python shell inside the MacBook terminal),
the final value of models has a equations_ of 9 equations (so 9 rows).
The model.equations_[8, :] (so the last equation found by the model "model") is:
complexity 17
loss 0.00174
score 1.872243
equation '(((1.4335294 / (x0 + 1.6759317)) + 0.6227007) / (0.21702437 - ((-1.312221 - (-0.2316349 / x0)) - x0)))
sympy_format (0.6227007 + 1.4335294/(x0 + 1.6759317))/(x0 +...
lambda_format PySRFunction(X=>(0.6227007 + 1.4335294/(x0 + 1...
Name: 8, dtype: object

However, when I run this exact equation using a Julia helper object ("jl"), I get a infinity value for the loss.
Specifically, I use the "jl" julia helper and I run:

jl.eval("""
# YOU CAN USE THIS CODE BLOCK TO TEST OUT THE VALUE OF THE CUSTOM OBJECTIVE FUNCTION

function eval_loss_real(tree, dataset, options)
    MAX_THRESHOLD_ALLOWED = 10000
    

    prediction, flag = eval_tree_array(tree, dataset, options)
    !flag && return L(Inf)
    
    num_points = length(prediction)
    
    for i in 1:num_points
        cur_expr_val = prediction[i]
        if (cur_expr_val > MAX_THRESHOLD_ALLOWED)
            println("prediction was greater than MAX_THRESHOLD_ALLOWED")
            return (Inf)
        end
        if (cur_expr_val < 0)
            println("prediction was less than 0")
            return (Inf)
        end
    end
    
    
    
    
    # BELOW SECTION IS TO CHECK FOR IF THE GUESSED FUNCTION TAKES ON ANY POSITIVE INFINITY VALUES
    # OR NEGATIVE INFINITY VALUES ON THE DOMAIN X=[0, 1].
    
    
    # NOTE: reshape([x], 1, 1) should give same dimensions
    # as reshape([x;], 1, 1).
    # f(x) = x -> 1/(eval_tree_array(tree, reshape([Float32(x)], 1, 1), options))
    f(x) = 1 / ((eval_tree_array(tree, reshape([Float32(x)], 1, 1), options))[1][1])
    
    # fp(x) = ForwardDiff.derivative(f,float(x))
    fp(x) = ( -1 / ((eval_tree_array(tree, reshape([Float32(x)], 1, 1), options))[1][1])^2 * (eval_diff_tree_array(tree, reshape([Float32(x)], 1, 1), options, 1))[2][1])
    
    global found_a_recip_root = true
    global recip_roots = -3
    try
        global recip_roots = find_zero((f, fp), 0.3, Roots.Newton(), maxevals=100)
        #println("recip_roots is", recip_roots)
    catch
        global found_a_recip_root = false
        # "Convergence Failed"
    end
    
    # If the function f(x) represented by the tree "tree"
    # has a value f(x) close to positive infinity or negative infinity
    # when x is in [0, 1], then return L(Inf).
    if found_a_recip_root
        if (recip_roots >= 0 && recip_roots <= 1)
            return (Inf)
        end
    end
    
    
    # BELOW SECTION IS TO CHECK FOR IF THE GUESSED FUNCTION TAKES ON ANY NEGATIVE VALUES OR ZERO VALUES
    # ON THE DOMAIN X=[0, 1].
    f(x) = ((eval_tree_array(tree, reshape([Float32(x)], 1, 1), options))[1][1])
    
    fp(x) = ((eval_diff_tree_array(tree, reshape([Float32(x)], 1, 1), options, 1))[2][1])
    
    global found_a_root = true
    global roots = -3
    try
        global roots = find_zero((f, fp), 0.3, Roots.Newton(), maxevals=100)
        #println("roots is", roots)
    catch
        global found_a_root = false
        # "Convergence Failed"
    end
    
    if found_a_root
        if (roots >= 0 && roots <= 1)
            println("found a root so returned L(Inf) and string rep of the tree was")
            println(string_tree(tree, options))
            return (Inf)
        end
    end
    
    
    
    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.01)
        global integral = my_tuple[1]
        global err = my_tuple[2]
    catch
        return (Inf)
    end
    
    
    
    norm_constant = integral
    if (integral <= 0)
        println("integral was <=0")
        return (Inf)
    end
    
    if (isinf(integral))
        println("integral was Inf")
        return (Inf)
    end
    
    actual_probs = (prediction) / norm_constant
    
    # length(actual_probs) equals length(prediction)
    
    # We really care about the EXPECTED VALUE OF L, where L = -1 * sum(log.(actual_probs)) / (length(prediction))
    # given that dataset.X is generated by the true probability distribution.
    # E(L) is really equal to the definite integral from x=0 to x=1 of ln(guessed_f(x)) * f_true(x) dx.
    # As num_points goes to +infinity, the variance of L becomes 0.
    
    final_loss = exp(-1 * sum(log.(actual_probs)) / (length(prediction)))
    #println("final_loss is ", final_loss)
    return final_loss
end

options = SymbolicRegression.Options(
    binary_operators=[+, *, /, -, ^],
    unary_operators=[cos, sin],
    enable_autodiff=true
)

num_points = 100

x = rand(Beta(3, 2), num_points)
x = reshape(x, 1, num_points)
x = Float32.(x)

x1 = Node(;feature=1)

#tree = -0.828202055*(x1*x1) + x1 - 0.14140837

# tree = x1 + x1
# tree = x1

# Final_loss was 0.67114264 (when I allowed it to not return Infinity, and used a norm_constant of 1).
#tree = -0.094332316551433*x1 + 1.7470982*x1/(1.67670228876569*x1 - 0.282760834105208)

# On model.equations_, the loss was 0.001739
tree = 0.21702437*(0.6227007 + 1.4335294/(x1 + 1.6759317))/(x1 + 1.52924537 - 0.2316349/x1)

tree = Node{Float32}(tree)

println(string_tree(tree, options))

my_loss = eval_loss_real(tree, x, options)
println("FINAL LOSS was ", my_loss)
""")

The output of running this jl.eval() was:

WARNING: Method definition (::Main.var"#f#52"{tree, options})(Any) in module Main at none:26 overwritten at none:69.
WARNING: Method definition (::Main.var"#fp#53"{tree, options})(Any) in module Main at none:29 overwritten at none:71.
WARNING: Method definition (::Main.var"#f#56"{tree, options})(Any) in module Main at none:169 overwritten at none:196.
WARNING: Method definition (::Main.var"#fp#57"{tree, options})(Any) in module Main at none:172 overwritten at none:198.
((0.21702437 * (0.6227007 + (1.4335294 / (x1 + 1.6759317)))) / ((x1 + 1.5292454) - (0.2316349 / x1)))
prediction was less than 0
FINAL LOSS was Inf

And it's not just this one time that I ran model.fit(X, y) that I noticed a difference between the model's reported loss of the model's last-found function and the usually-more correct outputted value of passing the model's last-found function into a Julia objective function manually and manually running the jl.eval() of that Julia code.

Another question: What would "WARNING: Method definition (::Main.var"#f#52"{tree, options})(Any) in module Main at none:26 overwritten at none:69." mean, what are its negative consequences, and how do I change my code to avoid this warning message?

[MilesCranmer]

You can get the number of rows from dataset.n. Alternatively using size(dataset.X, 2). You don’t want to use length because it will count all elements if there are multiple features.

[unary-code]

Yeah, in the model's objective function, I just changed length(prediction) in the formula for the variable "final_loss" to sizeof(prediction)[1]. After doing println(sizeof(prediction)[1]), I figured out that that is 100, so that's what I want.
However, even after making this change, the model's reported value of loss is clearly wrong and different from the returned output of the objective function when using the Julia helper object "jl" (when passing in a particular math function).

See my most recent post for this new issue: #375

[MilesCranmer]

With so much Julia customization, it might be easier to work on this directly from Julia? i.e., with SymbolicRegression.jl?