Why direct run of Julia code gives different output than loss value reported by model.equations_

[unary-code]

A direct run of the objective function by the Julia Helper gave this answer:

Matrix{Float64}
(1, 100)
Vector{Float64}
(100,)
my_loss0.9999999999999999

when given this code:

jl.eval("""

function process_eval_tree(tree, x, options, take_recip)
    #println("RUN EVAL x=", x)

    #x_float = Float32(x)
    #println("string_tree(tree, options) is", string_tree(tree, options))
    func_val, func_complete = (eval_tree_array(tree, reshape([x], 1, 1), options))
    #println("func_val is ", func_val)
    #println("func_complete is ", func_complete)
    
    if !func_complete
       return error("Evaluation of tree failed.")
    end
    
    if take_recip && func_val[1] == 0
        return Inf
    end
    
    if take_recip
        return (1/func_val[1])
    end
    return func_val[1]
end

function process_eval_diff_tree(tree, x, options, take_recip)
    #x_float = Float32(x)
    
    #println("RUN DERIV x=", x)
    
    deriv_val, gradient, deriv_complete = (eval_diff_tree_array(tree, reshape([x], 1, 1), options, 1))
    if !deriv_complete
       return error("Evaluation of tree's derivative failed.")
    end
    
    #println("typeof and size of gradient are")
    #println(typeof(gradient))
    #println(size(gradient))
    
    if isinf(gradient[1])
        return Inf
    end
    
    if take_recip
        recip_func_val = process_eval_tree(tree, x, options, take_recip)
        #println("recip func_val is", recip_func_val)
        #println("gradient[1] is", gradient[1])
        if isinf(recip_func_val)
            return Inf
        end
        return ( -1 * (recip_func_val * recip_func_val) * (gradient[1]))
    end
    
    return gradient[1]
end

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

    x_init_1 = Float64(0.99)
    x_init_2 = Float64(0.01)
    
    MAX_EVALS = 30
    
    tree = Node{Float64}(tree)
    dataset_X = Float64.(dataset.X)
    
    # size(dataset.X) is (100, 1)
    # println("size(x) is ", size(dataset.X))

    prediction, flag = eval_tree_array(tree, dataset_X, options)
    !flag && return L(Inf)
    
    # size(prediction) is (100,)
    #println("size(prediction) is ", size(prediction))
    
    # num_points is 100
    num_points = size(prediction)[1]
    
    # println("num_points is ", num_points)
    
    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([x], 1, 1), options))
    f_recip(x) = process_eval_tree(tree, x, options, true)

    fp_recip(x) = process_eval_diff_tree(tree, x, options, true)
   
    global found_a_recip_root = false
    global recip_roots = -3
    global recip_roots_2 = -4
    try
        #global recip_roots = find_zero((f_recip, fp_recip), Float64(0.3), Roots.Newton(), maxevals=10)
        global recip_roots = find_zero(f_recip, x_init_1, Order1(), maxevals=MAX_EVALS)
        global found_a_recip_root = true
    catch
        my_var = 1
    end
    
    try
        #global recip_roots_2 = find_zero((f_recip, fp_recip), Float64(0.01), Roots.Newton(), maxevals=10)
        global recip_roots_2 = find_zero(f_recip, x_init_2, Order1(), maxevals=MAX_EVALS)
        global found_a_recip_root = true
    catch
        my_var = 1
    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
        if (recip_roots_2 >= 0 && recip_roots_2 <= 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) = process_eval_tree(tree, x, options, false)
    
    fp(x) = process_eval_diff_tree(tree, x, options, false)
    
    global found_a_root = false
    global roots = -5
    global roots_2 = -6
    try
        #global roots = find_zero((f, fp), Float64(0.3), Roots.Newton(), maxevals=100)
        global roots = find_zero(f, x_init_1, Order1(), maxevals=MAX_EVALS)
        global found_a_root = true
        #println("roots is", roots)
    catch
        my_var = 1
        # "Convergence Failed"
    end

    try
        #global roots = find_zero((f, fp), Float64(0.3), Roots.Newton(), maxevals=100)
        global roots_2 = find_zero(f, x_init_2, Order1(), maxevals=MAX_EVALS)
        global found_a_root = true
        #println("roots is", roots)
    catch
        my_var = 1
        # "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
        if (roots_2 >= 0 && roots_2 <= 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 -> begin
       output, completion = (eval_tree_array(tree, reshape([x], 1, 1), options))
       if !completion
           error("Integration failed.")
       end
       output[1]
    end, Float64(0), Float64(1), rtol=0.1)
        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.
    
    # println("num_points is ", num_points)
    
    final_loss = exp(-1 * sum(log.(actual_probs)) / (num_points))
    if final_loss < 0.1
        println("string rep of tree is ", string_tree(tree, options))
        println("norm_constant is ", norm_constant)
        println("final_loss is ", final_loss)
    end
    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 = Float64.(X)
#X = Float32.(X)
println(typeof(X))
println(size(X))

y = 1000. * ones(num_points)
# y = Float64.(y)
#y = Float32.(y)
println(typeof(y))
println(size(y))


my_dataset = Dataset(X, y)

x1 = Node(;feature=1)

tree = ((x1 + x1) / (0.19358382 * x1))

my_loss = eval_loss(tree, my_dataset, options)
println("my_loss", my_loss)
""")

However, the exact same code (this time, I made sure to copy the entire objective function inside the variable named "objective", when run by model.fit(X, y), results in a model.equations_ that reports a loss value of 0.084261 for the function "((x0 + x0) / (0.19358382 * x0))
".

model.equations_ does not report the lowest-found loss value after calling eval_loss multiple times, right?

from pysr import PySRRegressor


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

using Roots
using ForwardDiff
using QuadGK

function process_eval_tree(tree, x, options, take_recip)
    #println("RUN EVAL x=", x)

    #x_float = Float32(x)
    #println("string_tree(tree, options) is", string_tree(tree, options))
    func_val, func_complete = (eval_tree_array(tree, reshape([x], 1, 1), options))
    #println("func_val is ", func_val)
    #println("func_complete is ", func_complete)
    
    if !func_complete
       return error("Evaluation of tree failed.")
    end
    
    if take_recip && func_val[1] == 0
        return Inf
    end
    
    if take_recip
        return (1/func_val[1])
    end
    return func_val[1]
end

function process_eval_diff_tree(tree, x, options, take_recip)
    #x_float = Float32(x)
    
    #println("RUN DERIV x=", x)
    
    deriv_val, gradient, deriv_complete = (eval_diff_tree_array(tree, reshape([x], 1, 1), options, 1))
    if !deriv_complete
       return error("Evaluation of tree's derivative failed.")
    end
    
    #println("typeof and size of gradient are")
    #println(typeof(gradient))
    #println(size(gradient))
    
    if isinf(gradient[1])
        return Inf
    end
    
    if take_recip
        recip_func_val = process_eval_tree(tree, x, options, take_recip)
        #println("recip func_val is", recip_func_val)
        #println("gradient[1] is", gradient[1])
        if isinf(recip_func_val)
            return Inf
        end
        return ( -1 * (recip_func_val * recip_func_val) * (gradient[1]))
    end
    
    return gradient[1]
end

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

    x_init_1 = Float64(0.99)
    x_init_2 = Float64(0.01)
    
    MAX_EVALS = 30
    
    tree = Node{Float64}(tree)
    dataset_X = Float64.(dataset.X)
    
    # size(dataset.X) is (100, 1)
    # println("size(x) is ", size(dataset.X))

    prediction, flag = eval_tree_array(tree, dataset_X, options)
    !flag && return L(Inf)
    
    # size(prediction) is (100,)
    #println("size(prediction) is ", size(prediction))
    
    # num_points is 100
    num_points = size(prediction)[1]
    
    # println("num_points is ", num_points)
    
    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([x], 1, 1), options))
    f_recip(x) = process_eval_tree(tree, x, options, true)

    fp_recip(x) = process_eval_diff_tree(tree, x, options, true)
   
    global found_a_recip_root = false
    global recip_roots = -3
    global recip_roots_2 = -4
    try
        #global recip_roots = find_zero((f_recip, fp_recip), Float64(0.3), Roots.Newton(), maxevals=10)
        global recip_roots = find_zero(f_recip, x_init_1, Order1(), maxevals=MAX_EVALS)
        global found_a_recip_root = true
    catch
        my_var = 1
    end
    
    try
        #global recip_roots_2 = find_zero((f_recip, fp_recip), Float64(0.01), Roots.Newton(), maxevals=10)
        global recip_roots_2 = find_zero(f_recip, x_init_2, Order1(), maxevals=MAX_EVALS)
        global found_a_recip_root = true
    catch
        my_var = 1
    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
        if (recip_roots_2 >= 0 && recip_roots_2 <= 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) = process_eval_tree(tree, x, options, false)
    
    fp(x) = process_eval_diff_tree(tree, x, options, false)
    
    global found_a_root = false
    global roots = -5
    global roots_2 = -6
    try
        #global roots = find_zero((f, fp), Float64(0.3), Roots.Newton(), maxevals=100)
        global roots = find_zero(f, x_init_1, Order1(), maxevals=MAX_EVALS)
        global found_a_root = true
        #println("roots is", roots)
    catch
        my_var = 1
        # "Convergence Failed"
    end

    try
        #global roots = find_zero((f, fp), Float64(0.3), Roots.Newton(), maxevals=100)
        global roots_2 = find_zero(f, x_init_2, Order1(), maxevals=MAX_EVALS)
        global found_a_root = true
        #println("roots is", roots)
    catch
        my_var = 1
        # "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
        if (roots_2 >= 0 && roots_2 <= 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 -> begin
       output, completion = (eval_tree_array(tree, reshape([x], 1, 1), options))
       if !completion
           error("Integration failed.")
       end
       output[1]
    end, Float64(0), Float64(1), rtol=0.1)
        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.
    
    # println("num_points is ", num_points)
    
    final_loss = exp(-1 * sum(log.(actual_probs)) / (num_points))
    if final_loss < 0.1
        println("string rep of tree is ", string_tree(tree, options))
        println("norm_constant is ", norm_constant)
        println("final_loss is ", final_loss)
    end
    return final_loss
end
"""


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)",
    progress=True
)

import numpy as np
num_points = 100
# (1, 100) versus (100, 1)
# I think the intended shape was (1, 100), right?
X = np.random.beta(3, 2, (num_points, 1))

print(X.shape)
#num_points = X.shape[0] # num_points is 100

y = 1000 * np.ones(num_points)

print(X)

model.fit(X, y)

print(model.equations_)

Why would this difference in output happen between the loss value printed by the Julia helper "jl" object and the loss value reported by model.equations_ ?

Feel free to use https://www.diffchecker.com/ so you can quickly understand the very small differences between the Julia helper code and the objective function used by the model "model".

[unary-code]

Another example:

model.equations_.iloc[5, :] is:

i= 5 th equation found
complexity          11
loss          0.260501
score         0.584707
Name: 5, dtype: object
(((x0 / 0.2732506) - x0) + (0.2732506 - (x0 - 0.11180633)))
PySRFunction(X=>1.65964429721289*x0 + 0.38505693)

Output of the Julia Helper object jl's code is:

Matrix{Float32}
(1, 100)
Vector{Float32}
(100,)
my_loss 0.87469167

where the Julia Helper object jl's code (using the Float32 function 1.65964429721289*x0 + 0.38505693) is:

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 = Float64.(X)
X = Float32.(X)
println(typeof(X))
println(size(X))

y = 1000. * ones(num_points)
# y = Float64.(y)
y = Float32.(y)
println(typeof(y))
println(size(y))


my_dataset = Dataset(X, y)

x1 = Node(;feature=1)

#tree = ((x1 + x1) / (0.19358382 * x1))

#tree = (-69.8075844683988*x1 - 0.6917331)/(-1*x1 - 0.93741227) + 0.021739371
tree = 1.65964429721289*x1 + 0.38505693
tree = Node{Float32}(tree)

my_loss = eval_loss(tree, my_dataset, options)
println("my_loss", my_loss)

Edit: My Julia Helper object also prints out a accurate loss value of:

Matrix{Float32}
(1, 100)
Vector{Float32}
(100,)
my_loss0.89749175

when I use

tree = (((x1 / 0.2732506) - x1) + (0.2732506 - (x1 - 0.11180633)))

instead of the simplified format

tree = 1.65964429721289*x1 + 0.38505693

[MilesCranmer]

Hi @unary-code,

Could you summarize the question/issue? I don't quite follow.

Thanks,
Miles

[unary-code]

I'm wondering why model.equations_ has a loss value of 0.084261 for the function "((x0 + x0) / (0.19358382 * x0))", while the Julia helper object "jl" calculates a loss function of 0.9999999999999999, given the same exact code for the function eval_loss.
A different dataset X is used between the model.equations_ and the Julia helper object.

I know that the expected value of the loss is 1, when the function is a constant ("((x0 + x0) / (0.19358382 * x0))" is a constant with respect to x0) so I expect the loss value to be 1. So the Julia helper object returned a correct loss value, while the model calculated and reported a loss value of 0.084261, which is really inaccurate and has a tiny almost-zero probability of happening.

The expected value of the loss is e to ( the integral from x=0 to x=1 of (- ln(f_hat(x) / norm_constant) * 12 * (x^2) * (1-x)) dx).
where norm_constant is the integral of f_hat(x) from x=0 to x=1.

See the below picture to understand what I'm trying to say.

[MilesCranmer]

Oh. It’s because you aren’t passing the loss function to the Julia version. So it’s just doing MSE.

[unary-code]

The Julia helper object is doing MSE for its loss value?

[unary-code]

Oh. It’s because you aren’t passing the loss function to the Julia version. So it’s just doing MSE.

Could you explain more? You mean to the Julia helper object "jl" ?

[MilesCranmer]

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

Missing the loss function

[unary-code]

Wait, I did do "options = SymbolicRegression.Options(
binary_operators=[+, *, /, -, ^],
unary_operators=[cos, sin],
enable_autodiff=true
)"
in my Julia code that the Julia helper object ran.

If you look above in the top post I made,
"
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)
" is part of the code that my Julia helper object ran.

How are you intended to create a Options object?
https://astroautomata.com/SymbolicRegression.jl/stable/api/#Options

[unary-code]

I'm pretty sure that neither the Julia helper object or the model is doing MSE for the loss

[unary-code]

Interesting:

I believe that when I ran:

jl.eval("""
    options = SymbolicRegression.Options(;
    binary_operators=[+, *, /, -, ^],
    unary_operators=[cos, sin],
    enable_autodiff=true
    )
    x1 = Node(;feature=1)

    tree = (((x1 - (x1 * x1)) / (((-0.37249994 + 0.40882915) / x1) - (0.0046278937 * x1))) + (0.0046278937 * 0.28049845))

    err = -3
    integral = 1
    err = -3
    try
        my_tuple = quadgk(x -> begin
       output, completion = (eval_tree_array(tree, reshape([x], 1, 1), options))
       if !completion
           return L(Inf)
       end
       output[1]
    end, Float64(0), Float64(1), rtol=0.001)
        integral = my_tuple[1]
         err = my_tuple[2]
    catch
        return L(Inf)
    end

    println("integral is", integral)
    println("err is", err)
    """)

I got this error:

┌ Warning: Assignment to `integral` in soft scope is ambiguous because a global variable by the same name exists: `integral` will be treated as a new local. Disambiguate by using `local integral` to suppress this warning or `global integral` to assign to the existing global variable.
└ @ none:22
┌ Warning: Assignment to `err` in soft scope is ambiguous because a global variable by the same name exists: `err` will be treated as a new local. Disambiguate by using `local err` to suppress this warning or `global err` to assign to the existing global variable.
└ @ none:23

and this was printed out:

integral is1
err is-3

However, I changed "global integral = my_tuple[1]" to "integral = my_tuple[1]" in my model's objective function eval_loss,
and now the model actually found a equation that numerically is close to the true function 12*x^2 (1-x).
So perhaps this was the issue?

My thoughts for why the model's loss value was inaccurate:
I know it is because the value of the variable "integral" was inaccurate.

Since the PySR package creates multiple threads which call eval_loss at the same time, the value of the variable "integral" that was being used to divide "predictions" by, was actually coming from a different call to eval_loss with a different parameter tree from a different integral.

When I did "global integral = my_tuple[1]", this actually changed the value of integral for a different call to eval_loss function for a different parameter "tree".

Question: Is there any way to make the call to model.fit(X, y) only create one thread instead of multiple threads? A.K.A. is there a way to make the call to model.fit(X, y) run completely sequentially, instead of parallel?

Question: I did "progress=False" when I constructed the model, but when I ran model.fit(X, y), it still showed the progress bar. Any reason this happens?