[BUG] min binary function raises errors during prediction

[gdelab]

Describe the bug
When I use "min" and "max" in my binary operators, fitting goes well but prediction raises an error.

Version (please include the following information):

• OS: ubuntu 20.04
• Julia version 1.8.0
• Python version 3.8
• PySR version 0.10.1
• Does the bug still appear with the latest version of PySR? Yes

Configuration

import numpy as np
from pysr import PySRRegressor

data_x = np.random.randint(-100, 100, (200, 2))
data_y = np.minimum(data_x[:,0], 8)

model = PySRRegressor(
    model_selection="best",  # Result is mix of simplicity+accuracy
    niterations=40,
    binary_operators=[
        "min", "max"
    ],
    unary_operators=[
    ],
    extra_sympy_mappings={
    },
    # ^ Define operator for SymPy as well
    loss="loss(x, y) = (x - y)^2",
    warm_start=False,
)
print('Fitting data')
model.fit(data_x, data_y)
print(model)
print(model.get_best())

print(model.predict(data_x))

Error message
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Happens in predict() call. The fit result is correct.

Additional context
I tried adding {'min': np.minimum}, or {'min': min}, or {'min': sympy.Min} to extra_sympy_mappings. None worked.

I found a workaround :

for row in data_x: 
    model.predict(row)

Not very satisfying though.

[gdelab]

Thanks for the amazing package by the way !

[MilesCranmer]

Hi @gdelab,

Thanks for posting. I didn't that min and max would work out of the box, so it is cool that they do!

I wonder if the reason these don't work at the prediction step is because min is already a defined function in python, so sympy just tries to interpret it as that, regardless of what you set for extra_sympy_mappings. One thing you could try is using a different name for the functions:

model = PySRRegressor(
    binary_operators=[
        "mmin(x, y) = min(x, y)",
        "mmax(x, y) = max(x, y)",
    ],
    extra_sympy_mappings={"mmin": sympy.Min, "mmax": sympy.Max}
)

Let me know if this works. Everything else looks fine though.

Cheers,
Miles

[MilesCranmer]

Actually, wait. I just looked and sympy seems to be doing it's job. If you look at model.sympy(i) where i is the row of the equation, you can see that it's correctly parsed to sympy!

I think the reason it doesn't work in the predict stage is because sympy's lambdify function does not have a numpy mapping for np.minimum yet. In other words, this is a sympy bug.

I'll think of a hack around this...

[MilesCranmer]

So I did something similar here: #150

You can hack this with the sympy.Piecewise function, which does correctly map to numpy equivalents. To use this:

model = PySRRegressor(
    binary_operators=[
        "min", "max"
    ],
    extra_sympy_mappings={
        "min": lambda x, y: sympy.Piecewise((x, x < y), (y, True)),
        "max": lambda x, y: sympy.Piecewise((y, x < y), (x, True)),
    }
)

This correctly uses model.predict!

[MilesCranmer]

Would you be interested in submitting a PR to add min and max to the default sympy mappings library? i.e., here:

PySR/pysr/sr.py

Lines 40 to 78 in dc1a4e7

	sympy_mappings = {
	"div": lambda x, y: x / y,
	"mult": lambda x, y: x * y,
	"sqrt_abs": lambda x: sympy.sqrt(abs(x)),
	"square": lambda x: x**2,
	"cube": lambda x: x**3,
	"plus": lambda x, y: x + y,
	"sub": lambda x, y: x - y,
	"neg": lambda x: -x,
	"pow": lambda x, y: sympy.Function("unimplemented_pow")(x, y),
	"pow_abs": lambda x, y: abs(x) ** y,
	"cos": sympy.cos,
	"sin": sympy.sin,
	"tan": sympy.tan,
	"cosh": sympy.cosh,
	"sinh": sympy.sinh,
	"tanh": sympy.tanh,
	"exp": sympy.exp,
	"acos": sympy.acos,
	"asin": sympy.asin,
	"atan": sympy.atan,
	"acosh": lambda x: sympy.acosh(abs(x) + 1),
	"acosh_abs": lambda x: sympy.acosh(abs(x) + 1),
	"asinh": sympy.asinh,
	"atanh": lambda x: sympy.atanh(sympy.Mod(x + 1, 2) - 1),
	"atanh_clip": lambda x: sympy.atanh(sympy.Mod(x + 1, 2) - 1),
	"abs": abs,
	"mod": sympy.Mod,
	"erf": sympy.erf,
	"erfc": sympy.erfc,
	"log_abs": lambda x: sympy.log(abs(x)),
	"log10_abs": lambda x: sympy.log(abs(x), 10),
	"log2_abs": lambda x: sympy.log(abs(x), 2),
	"log1p_abs": lambda x: sympy.log(abs(x) + 1),
	"floor": sympy.floor,
	"ceil": sympy.ceiling,
	"sign": sympy.sign,
	"gamma": sympy.gamma,
	}

[gdelab]

Thanks for the thourough answer. Yes, sure for the PR !

[tanweer-mahdi]

@MilesCranmer thanks for your amazing answer! What is the status of this PR? I looked into the current default sympy mappings but couldn't find them implemented!

[MilesCranmer]

Yeah I don’t think it has been submitted yet. You’re more than welcome to if you want!

[tanweer-mahdi]

@MilesCranmer I cannot seem to find that sympy_mappings dictionary in the master branch. Looks like the code blob you shared has not been merged either. Do you want me to create all of these mappings?