Replacing solve() with solveset()

[asmeurer]

This is a tracking issue for the progress forward for replacing solve with solveset. We need to figure out what this will look like in the long run, as well as what should be done before the 1.0 release.

Some things:

• [ ] Replacing internal calls with solveset https://github.com/sympy/sympy/issues/8711
• [x] Add solveset (and what other functions?) to __init__.py
• [ ] Is the solveset API stable? Do any changes need to be made before 1.0?
• [ ] Figure out how to deprecate solve in favor of solveset. A suggestion was to make solve a higher-level wrapper to solveset. Deprecation may be needed regardless, however, because of API changes.
• [x] Document the deprecation (even if we don't end up doing a formal deprecation by 1.0, we should still document that the change is happening) A couple of things that needs to be done to make it, at par with solve:

implement the following

• [ ] Non-linear multivariate system
• [ ] Equations solvable by LambertW (Transcendental equation solver).
• [ ] nested (trig) expression e.g. issue #10217

[aktech]

A couple of things that needs to be done to make it, at par with solve:

• Non-linear multivariate system
• Equations solvable by LambertW (Transcendental equation solver).

after these two things are implemented, we are ready to completely replace solve.

[asmeurer]

Which functions from solveset should be added to __init__.py?

[aktech]

Which functions from solveset should be added to init.py?

I would say just solveset(), linsolve() & linear_eq_to_matrix.

[asmeurer]

See https://github.com/sympy/sympy/pull/10049. Let me know if anything else should be added.

[asmeurer]

https://github.com/sympy/sympy/pull/10049 was merged.

@hargup @aktech would love to hear your thoughts on the above checkboxes. I don't think https://github.com/sympy/sympy/issues/8711 necessarily needs to be fixed before the release, unless we do end up deprecating solve (probably not going to happen, unless you have LambertW and non-linear multivariate ready to go).

So what really needs to be answered is:

• Are there any API changes that need to be made, especially to the now public functions solveset, linsolve, and linear_eq_to_matrix?
• How will we deprecate solve? What do you think of my idea of keeping solve as a higher-level wrapper to solveset (e.g., solve would check assumptions on the symbols and pass domain automatically)?

[asmeurer]

After a video chat with @hargup and @aktech, we've decided not to do any full deprecations of solve for this release but just document that it will be replaced and that people should use solveset. We also agreed to replace solve with solveset in places like the tutorial.

[asmeurer]

@hargup @aktech do either of you want to write this documentation, or should I do it?

[aktech]

@hargup @aktech do either of you want to write this documentation, or should I do it?

I will take care of it.

[aktech]

@asmeurer @hargup See #10131 Let me know if anything else should be added.

[asmeurer]

That pull request was merged. I think we decided not to do any deprecations yet. So I will remove the milestone from this issue. If there is anything that still needs to be done before 1.0 let me know.

[dnabanita7]

Hi,I am new to SymPy.I am interested in extending solveset project-idea.I am participating in GSoC-2019.So,any leads to start working will be helpful.Thanks in advance.

[Shekharrajak]

Hi @Naba7 ,

You can find the discussions. progress report and blogs in Solvers idea description. Go through it and if you have any question, feel free to ask/disucss in mailing thread or comment.

[dnabanita7]

Thanks @Shekharrajak

[Mohitbalwani26]

I am planning to participate in GSOC this year and I have gone through previous years progress of implementing solveset as almsot all type of equations are implemented already like transcendental equations,lambert type equations. so should i focus on improving them or should i find other things which could be added to solveset?

[oscarbenjamin]

We don't have a multivariate solveset. We have linsolve and nonlinsolve but neither of them returns a proper set in the case of underdetermined systems.

Also it might be worth figuring out how to use the output of solveset. There are lots of places where solve is used but solveset might potentially be better except that I don't really know what to do with the ouput of solveset. If I get, say, a union of imagesets how do I use that to do anything afterwards?

[oscarbenjamin]

Btw since @Mohitbalwani26 since your other PR was for the ode module it's worth mentioning that a lot of work is needed for systems of ODEs.

[moorepants]

Also it might be worth figuring out how to use the output of solveset. There are lots of places where solve is used but solveset might potentially be better except that I don't really know what to do with the ouput of solveset. If I get, say, a union of imagesets how do I use that to do anything afterwards?

I've attempted to use the new solvers over the years and always found that their outputs were not practical to work with. In general, sets are annoying to work with and my code that uses them seems to be full of converting them to lists. And then for the solve results, most uses cases are looking for more elementary results.

[hargup]

I've attempted to use the new solvers over the years and always found that their outputs were not practical to work with.

This was a huge design miss, mostly a result of inexperience on my part. If you see the original plan, it talks about creating this idealistic output type but it never asks the question "Will this be of any use to the end user?"

Before moving forward with adding any more new functionality we should focus on making sure existing solveset functionality can be used well.

[Mohitbalwani26]

@oscarbenjamin I dont have experience regardiny any PR related to solveset but I am willing to learn so should i do work on solveset or ode(in which I have experience as i have gone through the module thoroughly)

[asmeurer]

If solveset returns a FiniteSet, you can just iterate over it. Similarly, if it returns an infinite solution set, you can iterate the solutions. The issue is one, if an expression cannot be solved, it just returns an ImageSet representing the solutions symbolically, and two, if a solution depends on a condition, then iterating over the cases of that condition may not be easy (perhaps sets need to be better about iterating Piecewise in those cases). I think the worst cases of solveset are when it returns a nasty looking image set, but that really just means it couldn't find the solution. Or if it finds some solutions but can't be sure that it is all of them. So maybe we need some way of filtering those out.

[oscarbenjamin]

@Mohitbalwani26: You should work on what interests you or what you know more about or what you want to learn more about. Personally I'm more interested in the ODE module than solveset so that's my opinion :)

[oscarbenjamin]

@hargup thanks for the background. I think that solveset is good and can still be made more useful though :)

@asmeurer I guess what you mean is that if the equation can't be solved you get a ConditionSet. An ImageSet is fine but the expression from an ImageSet is probably more useful than the set itself e.g.:

In [1]: solveset(sin(x))
Out[1]: {2⋅n⋅π | n ∊ ℤ} ∪ {2⋅n⋅π + π | n ∊ ℤ}

Here it would often be more useful to get:

In [2]: solveset(sin(x))
Out[2]: π⋅n

I think that solveset is good in that it gives formally valid results. A ConditionSet represents a case that is indeterminate or does not have closed form solutions and that's fine because there will always need to be cases like that.

The real question is what do users want?

Sometimes I use solve just because I want to know if a solution exists. For purely numeric univariate equations solveset is quite good here and I can use solveset(eq).is_empty to get a tri-valued result. However when there are undefined symbols solveset suffers the same problems as solve (#16861).

More often I want to use the expressions returned by solve e.g. by substituting them into something. For that case solveset is suboptimal but I think it could be adapted. At least it would be possible to make a helper that extracts solution expressions from unions of finite sets and imagesets. Intersections are harder but again I guess that represents another possible use case: maybe I just want a complete list of what might be a solution rather than a precise minimal list of solutions.

[moorepants]

More often I want to use the expressions returned by solve e.g. by substituting them into something. For that case solveset is suboptimal but I think it could be adapted.

This is the most common use case for me and my students.

[oscarbenjamin]

More often I want to use the expressions returned by solve e.g. by substituting them into something. For that case solveset is suboptimal but I think it could be adapted.

This is the most common use case for me and my students.

The question is what do you want to know about the solutions that are returned. Is it important that get a complete solution set or do you just want any solution (or as many as possible...)?

A common use case for me is just wanting any solution, ideally the simplest. Consider for example trying to find the coefficients in the method of undetermined coefficients. For that sort of case I wrote a function that can get an element from a possibly infinite Set that can wrap solveset:

from sympy import *
from sympy.testing.pytest import raises

class Get1Exception(Exception):
    """Raised if solveset1 does not return a single solution"""
    pass

class IndeterminateSet(Get1Exception):
    """Raised if solveset is unable to solve the equation"""
    pass

class EmptySetError(Get1Exception):
    """Raised if solveset finds that there are no solutions"""
    pass

def best(sols):
    """Return the simplest element of iterable sols"""
    return min(sols, key=count_ops)

def get1(aset):
    """Get any element of a Set"""
    if aset.is_empty:
        raise EmptySetError
    elif aset.is_empty is None or isinstance(aset, ConditionSet):
        raise IndeterminateSet
    elif aset.contains(0) is S.true:
        return S.Zero
    elif aset.contains(1) is S.true:
        return S.One
    elif aset.contains(-1) is S.true:
        return S.NegativeOne
    elif isinstance(aset, FiniteSet):
        return best(aset)
    elif isinstance(aset, Union):
        sols = set()
        for s in aset.args:
            try:
                sols.add(get1(s))
            except Get1Exception:
                pass
        if not sols:
            raise NotFound
        return best(sols)
    elif isinstance(aset, ImageSet):
        return aset.lamda(*get1(aset.base_pset))
    elif isinstance(aset, ProductSet):
        return tuple(get1(s) for s in aset.sets)
    else:
        raise NotImplementedError  # Add more cases...

def solveset1(f, symbol=None, domain=Complexes):
    '''Get any solution of an equation'''
    return get1(solveset(f, symbol, domain))

raises(EmptySetError, lambda: solveset1(x**2 + 1, x, Reals))
raises(IndeterminateSet, lambda: solveset1(x*exp(x)+sin(x), x, Reals))

assert solveset1(x**2 - 1, x) == 1
assert solveset1(sin(x)) == 0
assert solveset1(sin(x)-S.Half) == pi/6
assert solveset1(sin(x)-1) == pi/2

I also needed this diff from master (adding is_empty to ImageSet):

diff --git a/sympy/sets/fancysets.py b/sympy/sets/fancysets.py
index c04c46243f..5b62ffdcfd 100644
--- a/sympy/sets/fancysets.py
+++ b/sympy/sets/fancysets.py
@@ -405,6 +405,10 @@ def __iter__(self):
     def _is_multivariate(self):
         return len(self.lamda.variables) > 1

+    @property
+    def is_empty(self):
+        return self.base_pset.is_empty
+
     def _contains(self, other):
         from sympy.solvers.solveset import _solveset_multi

With that you can do:

In [10]: solveset1(x**2+1, x, Complexes)                                                                                                       
Out[10]: ⅈ

In [11]: solveset1(x**2-1, x, Complexes)                                                                                                       
Out[11]: 1

In [12]: solveset1(x*exp(x), x, Complexes)                                                                                                     
Out[12]: 0

In [13]: solveset1(sin(x-2), x, Complexes).simplify()                                                                                          
Out[13]: 2 - π

[Shekharrajak]

If we want output type something like solve, then we have solvify.

Pasting the solvify docstring

Solves an equation using solveset and returns the solution in accordance
    with the `solve` output API.
    Returns
    =======
    We classify the output based on the type of solution returned by `solveset`.
    Solution    |    Output
    ----------------------------------------
    FiniteSet   | list
    ImageSet,   | list (if `f` is periodic)
    Union       |
    EmptySet    | empty list
    Others      | None
    Raises
    ======
    NotImplementedError
        A ConditionSet is the input.
    Examples
    ========
    >>> from sympy.solvers.solveset import solvify, solveset
    >>> from sympy.abc import x
    >>> from sympy import S, tan, sin, exp
    >>> solvify(x**2 - 9, x, S.Reals)
    [-3, 3]
    >>> solvify(sin(x) - 1, x, S.Reals)
    [pi/2]
    >>> solvify(tan(x), x, S.Reals)
    [0]
    >>> solvify(exp(x) - 1, x, S.Complexes)
    >>> solvify(exp(x) - 1, x, S.Reals)
    [0]

[gschintgen]

For my applications I definitely prefer solveset's general approach with its more complete and detailed output. But obviously it quickly becomes frustrating as soon as e.g. trigonometry becomes involved and one wants to do some further processing with the obtained solutions.

Solvify is very valuable in this respect, but from my perspective it has a design defect: it replaces solveset. So if I want to first look at the complete solution given by solveset (including e.g. conditions and periodicity) and then use those solutions, I'd have to call solvify with the same input equation and have Sympy do all of its processing all over again. (Even for apparently simple trig equation this can take multiple seconds on a weaker machine.)

I had a look at solvify and it does what I suspected: the very first line is just a call to solveset.

Considering that the output of solveset has been annoying users for years (don't get me wrong: its output is also one of its main assets!), why don't we just provide some kind of "solvify helper"? Like this maybe: 1) Move out all of solvify's code to a helper, maybe solvify_set or solveset_extract or extract_solutions (opinions?), and reduce solvify to a basic wrapper that just calls solveset and then the new helper. 2) Make this "helper" a fully public function and include it in Sympy's main __init__.py so that it's directly available in interactive sessions with from sympy import *. 3) Actively advertise this helper, in particular in the Tutorials with examples showing a convenient workflow:

[In] solveset(some trig equation)
[Out] some union of imagesets
[In] solvify_set(_)
[Out] [base solutions of the imageset expressions]

This new convenience function could possibly be extended with optional keywords. (dict for having dictionaries as output or periodic to have e.g. pi/3+n*2*pi instead of just pi/3.) But those would be details.

[oscarbenjamin]

Considering that the output of solveset has been annoying users for years (don't get me wrong: its output is also one of its main assets!), why don't we just provide some kind of "solvify helper"?

That seems reasonable to me.

This could be like Mathematica's Reduce and ToRules.

https://reference.wolfram.com/language/ref/Reduce.html https://reference.wolfram.com/language/ref/ToRules.html

Loosely you can think of the equivalence:

Solve[eq, x] == ToRules[Reduce[eq, x]]

This is not a precise equivalence though because Solve will not return solutions where the variables solved for are not x.

The primary difference as compared to solveset though is that Reduce turns a boolean statement into an equivalent boolean statement e.g.:

>>> reduce(x**2 - 1, x)
x = 1 | x = -1
>>> torules(_)
[{x: 1}, {x:-1}]

The reason for this separation is that reduce can maintain a clear validity guarantee: the output is always logically equivalent to the input. Turning that into substitution rules cannot preserve that guarantee e.g. there might be conditions on symbols other than the unknowns:

>>> reduce(a*x = b, x)
((a != 0) & (x = b/a)) | ((a = 0) & (b = 0))

We can't extract and use the substitution rule {x: b/a} from here without losing the information that this is only valid if a != 0.

Generally I think it would be better if solveset used boolean statements rather than sets but this also demonstrates another possible helper for the output of solveset: there could be a function for turning the sets into boolean statements. Boolean statements are better because they can represent many things much more naturally e.g. the output of reduce(a*x = b, x) above in which some of the statements are conditions that do not involve the unknown x. The reason that ToRules can take the output of Reduce directly is because it can extract the symbol names from the boolean statements.

This new convenience function could possibly be extended with optional keywords. (dict for having dictionaries as output or periodic to have e.g. pi/3+n*2*pi instead of just pi/3.) But those would be details.

I don't think that this is just a detail. This is actually a significant limitation in solvify: by not returning a full solution set it negates the possible advantage of using solveset over solve in the first place. Both solve and solvify are unusable in many situations:

In [5]: x = symbols('x', real=True)

In [6]: integrate(abs(sin(x)), x)
Out[6]: 
⎧  cos(x)    for x ≤ 0
⎪                     
⎨2 - cos(x)  for x ≤ π
⎪                     
⎩cos(x) + 4  otherwise

Here integrate needs to call solve to find out when sin(x) = 0 so that it can split the output into regions where sin(x) <= 0 and where sin(x) >= 0. It uses solve:

In [7]: solve(sin(x))
Out[7]: [0, π]

Now integrate thinks that there are only three regions of interest x < 0, 0 < x < pi and pi < x but this is incorrect. Using solveset with its complete representation of the solution set would be better but the output of solveset is difficult to use:

In [8]: solveset(sin(x))
Out[8]: {2⋅n⋅π │ n ∊ ℤ} ∪ {2⋅n⋅π + π │ n ∊ ℤ}

The solvify function aims to make it usable but in doing so brings back precisely the same problem of using solve in the first place:

In [12]: from sympy.solvers.solveset import solvify

In [13]: solvify(sin(x), x, Reals)
Out[13]: [0, π]

What is the point of replacing a call to solve with a call to solvify if by design solvify has the same fundamental flaws as solve?

It might not be convenient for interactive use but for internal/programmatic use possibly the most useful way to represent the solution sets is as a union of sets where each is both an ImageSet and a ConditionSet:

solve(sin(a*x), x)
[
    # x = n*pi for n in Integers if a != 0:
    ({x: n*pi}, {n: Integers}, a != 0),
    # x = anything if a = 0:
    ({}, {}, a = 0),
]

This is analogous to what is returned whenusing Matlab's ReturnConditions parameter for its solve function: https://uk.mathworks.com/help/symbolic/solve-an-algebraic-equation.html

[gschintgen]

possibly the most useful way to represent the solution sets is as a union of sets where each is both an ImageSet and a ConditionSet

I agree that for programmatic use it'd be important to have a consistent output including all of the conditions. I'm finding your proposed output intriguing so I put together a proof of concept.

from sympy import *
from sympy.logic.boolalg import simplify_logic

def newsolvify(f, symbol, domain, solveset_out=None):
    """Solveset with alternative output. Proof of concept.

    https://github.com/sympy/sympy/issues/10006#issuecomment-1843359121
    """
    if not solveset_out:
        solveset_out = solveset(f, symbol, domain)

    if solveset_out is S.EmptySet:
        return []

    def solset_to_list(solset):
        soltuples = []
        if isinstance(solset, Union):
            for us in solset.args:
                soltuples += solset_to_list(us)

        if isinstance(solset, FiniteSet):
            soltuples += [ ({symbol: s}, {}, True) for s in solset ]

        elif isinstance(solset, ImageSet):
            baseset = solset.base_set
            expr = solset.lamda.expr
            assert len(solset.lamda.variables) == 1
            var = solset.lamda.variables[0]
            soltuples += [ ({symbol: expr}, {var: baseset}, True) ]

        elif isinstance(solset, ConditionSet):
            baseset = solset.base_set
            var = solset.sym
            cond = solset.condition
            if isinstance(baseset, FiniteSet):
                csym = list(cond.free_symbols)[0]
                cond = cond.as_set().as_relational(csym)
                soltuples += [ ({symbol: e}, {}, cond) for e in baseset ]
            elif baseset == S.Reals:
                soltuples += [ ({}, {}, cond)]
            else:
                y = Dummy('y')
                soltuples += [ ({symbol: y}, {y: baseset}, cond) ]

        elif isinstance(solset, Complement):
            baseset = solset.args[0]
            condset = solset.args[1]
            if isinstance(condset, ConditionSet):
                if condset.base_set == S.Reals:
                    cond = ~condset.condition
            else:
                cond = ~ (Contains(symbol, condset).as_set().as_relational(symbol))
                cond = simplify_logic(cond)
            basetuples = solset_to_list(baseset)
            soltuples += [ (bt[0], bt[1], bt[2] & cond) for bt in basetuples ]   

        return soltuples

    return solset_to_list(solveset_out)

Here are some of the inputs that I used to debug it:

newsolvify((x-1)*(x-2)*(x-3)*sin(3*x), x, S.Reals)
newsolvify((x-1)*(x-2)*(x-3)*sin(3*x)/(x**2-4), x, S.Reals)
newsolvify(x+exp(x), x, S.Reals)
newsolvify(abs(x)-a, x, Reals)
newsolvify((x-1)*(x-2)*sin(3*x)/(x-exp(x)), x, S.Reals)

That last input returns:

[({x: 1}, {}, True),
 ({x: 2}, {}, True),
 ({x: 2*_n*pi/3}, {_n: Integers}, Ne(x - exp(x), 0)),
 ({x: 2*_n*pi/3 + pi/3}, {_n: Integers}, Ne(x - exp(x), 0))]

As for interactive use, I think I'd prefer a simpler solvify type wrapper. I also realized that solvify can't just be modified to take solveset's output as sole input, since that would lose the variable name that is necessary to construct a meaningful subs dictionary.

[oscarbenjamin]

As for interactive use, I think I'd prefer a simpler solvify type wrapper.

I agree but it is better to start by having a completely correct general format and then it is easy to build something more like solvify on top. That is sort of the idea of solveset/solvify already but:

• solveset is not completely correct because it does not handle degenerate cases like a=0 in a*x = b.
• The output of solveset needs to be converted into a more structured format before it can be used for anything (it is already not a convenient format for solvify).
• Some of the choices made by solvify lead to the same problems as solve in the first place again.

None of these is really the biggest issue with the idea of replacing solve with solveset which is that solveset is only for univariate equations. It is a mistake to design a solver that is only for univariate equations because even univariate equations are often best handled by introducing new variables and turning them into systems of equations.