Bug: Solveset with piecewise constant functions

[trenta3]

Hi! I couldn't find a procedure to call for inverting a specified real-valued function.
For inverting I mean that it should return the different branches of the inverse, i.e. if I would ask it to invert $f(x) = x^2$ I would like to have returned $(g_1(y) = \sqrt{y}, y \in [0, \infty])$ and $(g_2(y) = - \sqrt{y}, y \in [0, \infty])$.

The "trick" of using solveset as in

from sympy import *
x, y = var("x y", real=True)
solveset(y - f(x), x, Reals)

seems not to specify the validity domain of the returned functions (i.e. in the example above does not return the result as Image(x -> sqrt(x), [0, \infty]) but only returns sqrt(x), sqrt(-x)).

Moreover it also appears to have major difficulties in dealing with the inversion of piecewise constant functions, such as the Heaviside function:

H = Lambda(x, Max(x, 0).diff(x))
solveset(y - H(x), x, Reals)

returns emptyset, while I would expect it to return a piecewise definition such as:

• If y = 0, then (-\infty, 0]
• If y = 1, then (0, \infty)
• Otherwise, Emptyset

Is there any function that does this that I have overlooked? How difficult would it be to integrate such functionality (returning more appropriate information from solveset in functional form) into sympy?

[MaqOwais]

For invert see (https://github.com/sympy/sympy/blob/70381f282f2d9d039da860e391fe51649df2779d/sympy/solvers/solveset.py#L176-L181
(documentation part)

I think for more vivid solution , it's under development :)

[trenta3]

Thanks a lot for the reference to sympy.solvers.solveset.invert_real.
I think the best thing is that it does not return incorrect result with the Heaviside function, despite being not able to simplify it.

[asmeurer]

You can also use plain solve. For Heaviside, I think it isn't implemented yet.

[trenta3]

Actually, looking into the code I found the source of the error, so I would like to "promote" this issue to a bug report.

>>> solveset(y - Heaviside(x), x, Reals)
EmptySet

I manually executed the lines of the function solveset in the current main branch.
After skipping the first checks, the lines that are executed are:

sympy/sympy/solvers/solveset.py

Lines 2026 to 2034 in 0e22875

	f, mask = _masked(f, Abs)
	f = f.rewrite(Piecewise) # everything that's not an Abs
	for d, e in mask:
	# everything *in* an Abs
	e = e.func(e.args[0].rewrite(Piecewise))
	f = f.xreplace({d: e})
	f = piecewise_fold(f)
	return _solveset(f, symbol, domain, _check=True)

in which y - Heaviside(x) is rewritten as Piecewise((y, x < 0), (y - Heaviside(0), Eq(x, 0)), (y - 1, x > 0)) and then is passed to _solveset.
It first skips all checks until this one:

sympy/sympy/solvers/solveset.py

Lines 934 to 941 in 0e22875

	elif f.is_Piecewise:
	result = EmptySet
	expr_set_pairs = f.as_expr_set_pairs(domain)
	for (expr, in_set) in expr_set_pairs:
	if in_set.is_Relational:
	in_set = in_set.as_set()
	solns = solver(expr, symbol, in_set)
	result += solns

in which it is split into pieces, and thus

solver(expr=y, symbol=x, in_set=Interval.open(-oo, 0)) 

is called with the first piece of the piecewise function, where solver = _solveset.
This time the code recognizes that expr y does not contain the symbol x and therefore outputs EmptySet.

sympy/sympy/solvers/solveset.py

Lines 915 to 918 in 0e22875

	if f.expand().is_zero:
	return domain
	elif not f.has(symbol):
	return EmptySet

Other pieces follow the same faith.

I think that one would like to return a ConditionSet(x = y) in this case, which at least would give meaningful results in this case, but it might break lots of stuffs elsewhere.
Another proposal would be to add a silent_variables list to _solveset to remember on which variables the function does depend, but is constant with respect to them.

I would like to submit a Pull Request for this, since I need it such functionality, and also because I would like to familiarize myself with sympy repo.
Let me know if I can do it and exactly how I could fix it.

[asmeurer]

You're of course free to submit a pull request. I don't know what the fix should look like. Someone more familiar with the solvers code will have to comment.

[trenta3]

Yes, I will wait for someone who known the solvers code to answer.
In the meanwhile I substituted the return EmptySet line in

sympy/sympy/solvers/solveset.py

Lines 915 to 918 in 0e22875

	if f.expand().is_zero:
	return domain
	elif not f.has(symbol):
	return EmptySet

with a return ConditionSet(symbol, Eq(f, 0), domain) and no test broke.

It should now only be a matter of deciding how to return such results as sets, since they involve more than one variable.

[gschintgen]

Yes, I will wait for someone who known the solvers code to answer.
In the meanwhile I substituted the return EmptySet line in

sympy/sympy/solvers/solveset.py

Lines 915 to 918 in 0e22875

	if f.expand().is_zero:
	return domain
	elif not f.has(symbol):
	return EmptySet

with a return ConditionSet(symbol, Eq(f, 0), domain) and no test broke.

It should now only be a matter of deciding how to return such results as sets, since they involve more than one variable.

This reminds me of #16861 and in particular #16861 (comment).