Longer runtime when running several PySR runs after another

[BrotherHa]

Hello,

I am using PySR for my master's thesis and I came across some strange behavior, that I have not found any comments about on this Github page yet.
When I run several PySR runs after another in one Skript (for example in a for-loop), the runs tend to take longer and longer, even though their settings are identical, no warmstart is enabled and I am creating a new PySRRegressor instance for each run.
An example function, that causes the behavior is given below together with the recorded increase in runtime and RAM Usage. The RAM usage was my first guess, but it is not increasing for the searches (see record). Also the clock frequency of the CPU is staying constant.
In the PySR output it is noticeable, that the number of evaluations per second is decreasing.
The increase in runtime is becoming higher when running a more extensive regression for example with many populations.
I really do not now, what might cause this issue and thereby I would be very grateful for your support!
Thank you very much in advance already.

Best regards

---

import numpy as np
import pandas as pd
import time
import tracemalloc
from pysr import PySRRegressor


# Training data and loss_function are created before

# PySR Settings
binary_operators = ["+", "-", "*", "/", "^"]
unary_operators = ["exp", "log", "sin", "atan", "sqrt", "square", "neg", "abs"]
complexity_of_operators = {"exp":2, "log":2, "sin":2, "atan":2, "sqrt":2, "square":2, "abs":3}
constraints = {'^':(-1,2),'sin':(2),'atan':(2)}
nested_constraints = {'atan':{'atan':0,'sin':1},'sin':{'sin':0,'atan':1},'exp':{'log':0},'log':{'exp':0},'sqrt':{'square':0,'neg':0},'square':{'sqrt':0,'neg':0,'abs':0},'neg':{'neg':0},'abs':{'abs':0,'neg':0}}
batching = True
batch_size = 500

# Start multiple PySR runs
memory = []
t_run_s=[]
for i in range(6):
    # Create regression object
    regression_object = PySRRegressor(loss_function=loss_function, niterations=100, warmup_maxsize_by=0.5,
                                      population_size=100, populations=25, tournament_selection_n=5,
                                      model_selection='accuracy', binary_operators=binary_operators,
                                      unary_operators=unary_operators, constraints=constraints,
                                      nested_constraints=nested_constraints, batching=batching, batch_size=batch_size,
                                      bumper=True, turbo=True)
    t_start = time.time()
    tracemalloc.start()
    regression_object.fit(training_data[['x1','x2','x3']],training_data[['y']])
    memory.append(tracemalloc.get_traced_memory())
    tracemalloc.stop()
    t_finish = time.time()
    t_run_s.append(t_finish - t_start)
# Print RAM usage and duration
print('Memory (first value: current RAM usage (uninteresting), second: peak RAM usage in run (interesting)):')
print(memory)
print('Time: ')
print(t_run_s)


# Resulting output:
#Memory (first value: current RAM usage (uninteresting), second: peak RAM usage in run (interesting)):
#[(1224377, 2029137), (947962, 1740351), (1160301, 2339805), (2075128, 3449471), (1254039, 2223228), (3452654, 3933324)]
#Time:
#[470.570111989975, 513.3470871448517, 609.9849309921265, 719.835132598877, 870.5115010738373, 1025.4808928966522]

[MilesCranmer]

Maybe related to #764? This turned out to be a Julia bug which I reported in JuliaLang/julia#56759 and is now fixed. Now we just wait for 1.11.3, or you can force it to use 1.10 (see PySR issue for how).

[BrotherHa]

Thank you very much for the quick answer! I tried using 1.10(.7), but unfortunately it did not solve the issue.

[MilesCranmer]

What is your loss_function?

[BrotherHa]

Thanks for your answer and your support! I am really not finding the reason and the whole issue seems quite weird to me. It also occurs with the PySR version 0.19.4, that is what I tried today. And the increase in run time gets bigger continuously for more runs and does not converge as far as I tried. I just came across this issue now, as I was not running many iterations before and the runs where not that extensive (e.g. smaller populations) and there the effect was far smaller.

My loss function is defined as follows. It depends on the particular input data and the options in the beginning are to identify the different y-variables, that all need differently shaped loss functions:

function custom_loss_function(tree, dataset::Dataset{T,L}, options, idx)::L where {T,L}
    if isapprox(dataset.y[1],0.0)
        p5 = 0.00813762
        p80 = 17.448
        a_1 = 7.838560195225847
        b_1 = 0.012206430080501547
        a_2 = 2
        b_3 = 6.53098971058287
        a_3 = 47.9579776388101
        c_3 = -117.475032351869
    elseif isapprox(dataset.y[1],0.05725155)
        p5 = 1.2804586
        p80 = 15.113608
        a_1 = 0.6248880798745418
        b_1 = 1.9206878219999999
        a_2 = 2
        b_3 = -6.09607886912224
        a_3 = 18.0350585617555
        c_3 = -9.43495030694149
    elseif isapprox(dataset.y[1],0.00017221528)
        p5 = 0.0040257196
        p80 = 1.8914138
        a_1 = 11.144567754841631
        b_1 = 0.006038579739964286
        a_2 = 2
        b_3 = 0.0400789096218153
        a_3 = 3.86298538267220
        c_3 = 1.23985086650155
    elseif isapprox(dataset.y[1],0.0063780397)
        p5 = 0.46299943
        p80 = 13.051656
        a_1 = 1.0391895552455737
        b_1 = 0.694499148
        a_2 = 2
        b_3 = 4.05385061247809
        a_3 = 34.2110131049562
        c_3 = -71.0354535190814
    end
    if idx == nothing
        y_pred, complete = eval_tree_array(tree, dataset.X, options)
        y = dataset.y
    else
        y_pred, complete = eval_tree_array(tree, dataset.X[:,idx], options)
        y = dataset.y[idx]
    end
    y_border = zeros(length(y))
    y_border[abs.(y).<p5] .= (a_1 * y[abs.(y).<p5]).^2 .+ b_1
    y_border[(abs.(y).>=p5).&(abs.(y).<=p80)] .= a_2 * abs.(y[(abs.(y).>=p5).&(abs.(y).<=p80)])
    y_border[abs.(y).>p80] .= a_3 * log.(abs.(y[abs.(y).>p80]) .+ b_3) .+ c_3
    if !complete
        return L(Inf)
    end
    if dataset.weights isa Array{<:Number,1}
        if idx == nothing
          custom_error = (sum( ( (-3 .* ( ( ( ((y_pred .- y) ./ (sqrt(2) .* y_border)).^2 ) .+ 1).^(-1) ) ) .+ 3) .* dataset.weights))/(sum(dataset.weights))
        else
          custom_error = (sum( ( (-3 .* ( ( ( ((y_pred .- y) ./ (sqrt(2) .* y_border)).^2 ) .+ 1).^(-1) ) ) .+ 3) .* dataset.weights[idx]))/(sum(dataset.weights[idx]))
        end
    else
        custom_error = (sum( (-3 .* ( ( ( ((y_pred .- y) ./ (sqrt(2) .* y_border)).^2 ) .+ 1).^(-1) ) ) .+ 3))/(length(y))
    end
    return custom_error
end

[MilesCranmer]

Very weird. If you do a time.sleep(100000) after the first search is complete, what does htop show? Is Julia still running? I wonder if it’s some zombie process still doing the search that wasn’t killed or something? Just trying to think of explanations.

How are you letting each run quit? Do they run until the end or do you quit them early?

Do you see this with bumper set to False?

[MilesCranmer]

Also, if you’re up for getting your hands dirty, you could try building Julia at this commit: JuliaLang/julia#56801, and then forcing PySR to use that Julia. (via juliapkg.require_julia before importing PySR). My gut feeling is that it’s due to the bug that was fixed in that Julia PR (which was also found on 1.10 apparently).

I’m curious if turning off bumper helps too though.

[BrotherHa]

Thanks for your answers! During the sleep, the python.exe of the process showed no activity in regards to memory or CPU at all. And there was also no process. that seemed suspicious to me. I checked that with ntop and the Process Explorer, as I am working on a Windows machine.
Each run quits when niterations is met. I do not quit them early.
And even with bumper set to False the behavior occurs.

Concerning your second answer I tried to run Julia at the state of the commit, but unfortunately I did not manage. How can I input the SHA of the commit into the require_julia function or is there another way?

What I do for now is to "outsource" the for loops I run by having a .bat-file, that calls a python-script to run PySR for several times and I save the parameters of each regression as .pkl and read it in the next iteration to continue where I ended. This works without increasing runtime for many iterations.

[BrotherHa]

I am running some iterations now on a server with many core and multiple CPUs and when turning parallelism to "multiprocessing" the behavior is not occurring any more. I really have no idea, what might be the reason for the issue, when running on a single machine. If you could specify how to force PySR to use a certain Julia commit I would still try to see, if that solves the problem.
Thank you very much for your support!

[MilesCranmer]

Here's how to do it by installing the nightly version of Julia.

First, install Juliaup:

curl -fsSL https://install.julialang.org | sh
# on windows: winget install julia -s msstore

Then, install the "nightly" channel:

juliaup add nightly

This is the most recent version of Julia in the master branch. (Note there is a chance it is broken, since it's not as stable as the released versions).

Then, when you start Python, set the PYTHON_JULIAPKG_EXE environment variable to point to the nightly Julia which you can find with this:

julia +nightly --startup-file=no -e 'println(Base.julia_cmd())'
# For me, I get: `/Users/mcranmer/.julia/juliaup/julia-nightly/bin/julia -C native -J/Users/mcranmer/.julia/juliaup/julia-nightly/lib/julia/sys.dylib -g1 --startup-file=no`

Set this when you start python:

import os
os.environ["PYTHON_JULIAPKG_EXE"] = "/Users/mcranmer/.julia/juliaup/julia-nightly/bin/julia"

# Then, import pysr
import pysr

And it will automatically install under this specific Julia version!