New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Rewriting the Univariate Solver #7523

Merged
merged 3 commits into from Dec 16, 2014

Conversation

Projects
None yet
8 participants
@hargup
Copy link
Member

hargup commented May 25, 2014

This PR aims to rewrite a cleaner and robust univariate equation solver.

  • polynomial
  • rational
  • real trigonometric
  • write 4 tests to complete coverage
  • correct typo
  • correct this failure

Nonblocking todo

  • functions solvable by LambertW
  • functions that can be recast as polynomials with a change of variables this, for example; this can be factored out of solve where multiple generators are handled
  • use something like this to handle the XFAILed test test_real_imag_splitting1, this will be handled in the set module.
@hargup

This comment has been minimized.

Copy link
Member

hargup commented May 25, 2014

@skirpichev @mrocklin Please have a look.

"""
Robust univariate equation solver
"""
if f.is_Mul:

This comment has been minimized.

@mrocklin

mrocklin May 25, 2014

Member

Is this relationship more general than for just univariate systems?

This comment has been minimized.

@hargup

hargup May 25, 2014

Member

Yes, I can't think of a reason why it won't work for a multivariate system. f(x, y)*g(x,y) = 0 implies either f(x, y) = 0 or g(x, y) = 0 or both.

for m in f.args:
solns = solve_univariate(m, symbol)
result.update(set(solns))
return list(result)

This comment has been minimized.

@mrocklin

mrocklin May 25, 2014

Member

Maybe this can be tightened up into a single comprehension? I'm not sure which is better for readability though.



def invert(f, symbol):
""" Returns the list of the inverse function """

This comment has been minimized.

@mrocklin

mrocklin May 25, 2014

Member

A simple doctest would be good here to help explain what this does.

# also ploting for the piecewise functions doesn't work, it wll be easy
# to implement create a issue for it.
# TODO: write the docstring for the as_expr_set_pairs method for
# piecewise functions

This comment has been minimized.

@mrocklin

mrocklin May 25, 2014

Member

I find that unindented comments interrupt the nested structure of the code. For me personally this makes things harder to read. I'm not sure what standard style is though.

"""
Solves the equation as polynomial or converting
it to a polynomial.
"""

This comment has been minimized.

@mrocklin

mrocklin May 25, 2014

Member

Again, I'm a big fan of doctests. It's a great introduction to the code.

Solves the equation as polynomial or converting
it to a polynomial.
"""
def solve_as_poly_gen_is_pow(poly):

This comment has been minimized.

@mrocklin

mrocklin May 25, 2014

Member

Can this function live outside of solve_as_poly? This will make it easier to test and also help to keep functions small.


if f.is_polynomial(symbol):
solns = roots(f, symbol)
no_roots = sum([solns[key] for key in solns.keys()])

This comment has been minimized.

@mrocklin

mrocklin May 25, 2014

Member

sum(solns.values()) ?



def test_invert():
assert invert(x + 3, x) == [x - 3]

This comment has been minimized.

@mrocklin

mrocklin May 25, 2014

Member

It might not be important, but I feel the need for an answer like y - 3, not x - 3. Is it strange to use the independent variable from the input as the independent variable in the output?

This comment has been minimized.

@hargup

hargup May 25, 2014

Member

If chose a different output variable from the input variable we will have to explicitly as for it as a parameter. The other option is dummy, but they look terrible in printing and they will also mess up the tests, I suppose.

This comment has been minimized.

@asmeurer

asmeurer May 25, 2014

Member

You can just allow the dummy to be passed in as an argument, like

def invert(eq, x, y=None):
    y = y or Dummy('y')

This comment has been minimized.

@smichr

smichr May 26, 2014

Member

From what I've seen if f(x) = 2x then f^-1(x)=x/2. x is just a variable and it is the function markup that indicates what will be returned for a given x. Here, your "markup" is the function name "invert". So I think using x is fine.

This comment has been minimized.

@asmeurer

asmeurer May 26, 2014

Member

My experience from a pedagogy point of view is that using the same symbols for a function and it's inverse is very confusing. Assumedly that confusion would carry over here as well.


assert invert(log(x), x) == [exp(x)]
assert invert(log(3*x), x) == [exp(x)/3]
assert invert(log(x + 3), x) == [exp(x) - 3]

This comment has been minimized.

@mrocklin

mrocklin May 25, 2014

Member

To motivate the return type of list there should probably be a case with multiple outputs?


@XFAIL
def test_fail_invert():
assert invert(x*log(x), x)

This comment has been minimized.

@mrocklin

mrocklin May 25, 2014

Member

What should this equal when it works?

return list(solns.keys())
else:
raise NotImplementedError
elif not f.is_Function and not f.is_Mul:

This comment has been minimized.

@mrocklin

mrocklin May 25, 2014

Member

Even though is_Foo are prevalent in SymPy I believe that current policy is to prefer isinstance when possible.

def solve_as_poly(f, symbol):
"""
Solves the equation as polynomial or converting
it to a polynomial.

This comment has been minimized.

@mrocklin

mrocklin May 25, 2014

Member

I don't understand this perfectly. Is the simple tagline "Solves a polynomial equation" correct? If not why not?

else:
raise NotImplementedError
else:
return solve_as_poly(f, symbol)

This comment has been minimized.

@mrocklin

mrocklin May 25, 2014

Member

Should there be a check here on f? For example consider the input solve_univariate(cos(x**2) - 1, x). This correctly raises a NotImplementeError but it does so from within solve_as_poly. Probably it should err sooner.

@mrocklin

This comment has been minimized.

Copy link
Member

mrocklin commented May 25, 2014

This looks good. Definitely easier for me to understand than old solve.

Two points of feedback:

  • I really value doctests, even just one simple one for each function that clearly outlines scope. Something as simple as the following is great
def solve_poly(f, symbol):
    """ Solve polynomial equations

    >>> solve_poly(x**2 - 1, x)
    [-1, 1]
    """
  • So far most of the logic is type checking. This suggests that we might try some _eval_solve method. I wouldn't do this yet of course (probably this is too restrictive), just thought I'd bring it up. Even if we don't do this it might still make sense to break up solve_univariate into smaller functions. We're not getting the reuse benefit of these functions (they're only likely to be called in solve_univariate, but we do get the benefit of breaking up testing into smaller units and documenting smaller chunks of code. My aesthetics are probably more extreme here than average (so you probably shouldn't listen to me), but if this were my project I would probably make solve_univariate look more like this:
def solve_univariate(expr, symbol):
    """
    ....
    """
    if f.is_polynomial(symbol):
        return solve_univariate_polynomial(expr, symbol
    elif not f.is_Function and not f.is_Mul:
        return solve_...(expr, symbol)

It has the benefit that you can read the high-level policy/cases easily and then grep for the parts that you care about. It's more like a hyper-linked webpage than a linearly laid out book. It's probably also easier to refactor when we change things in the future.

@mrocklin

This comment has been minimized.

Copy link
Member

mrocklin commented May 25, 2014

Are you rewriting solve logic or are you able to salvage code from current solvers.py?

It would be nice to satisfy a slightly modified version of the original test suite. You might consider copying over the old tests, removing the ones that don't make sense, XFAILing the others, and then slowly un-XFAILing them as they start to pass.

@hargup

This comment has been minimized.

Copy link
Member

hargup commented May 25, 2014

Are you rewriting solve logic or are you able to salvage code from current solvers.py?

I'm using the tests from the current solvers.py, the new solvers should be at least as good as the old one. I'm adding new test too, where ever I feel the current tests are not sufficient.

@hargup

This comment has been minimized.

Copy link
Member

hargup commented May 25, 2014

I don't understand this exactly. Am I right in assuming the following behavior?

>>> invert(Abs(x), x)
[-x, x]

Yes, that's the expected behavior of the function.

Is invert a cheap solve?

Yes, in some sense.

How far can we take this simple approach?

I don't think too far. I intend to keep it for simple case only and independent of the solve. This is intended for internal use only. Inverting a function is same as solving f(x) - y = 0, so solve will serve as general invert function.

If so then what about

>>> invert(cos(x), x)
?

Maybe this should return a set?

The same thing we should return for solve(cos(x) - a, x). Only sets seem to be a good candidate. btw I'm not incorporating trig functions in the initial re-factoring. Maybe after that I can change the output to set and after that we add code for trig functions.

Robust univariate equation solver
"""
if f.is_Mul:
result = set()

This comment has been minimized.

@skirpichev

skirpichev May 25, 2014

Contributor

move this to the beginning of the func?

@hargup hargup added GSoC labels May 25, 2014


def solve_univariate(f, symbol):
"""
Robust univariate equation solver

This comment has been minimized.

@skirpichev

skirpichev May 25, 2014

Contributor

just "univariate equation solver"

# We might dispach it into the functions themselves
if f.is_Symbol:
return [f]
if isinstance(f, exp):

This comment has been minimized.

@skirpichev

skirpichev May 25, 2014

Contributor

I don't like this growing if hell. I think, we have inverse methods in some function classes just for this purpose.

@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented May 25, 2014

I never thought _eval_solve made much sense as most methods would just be trivial, except for Add._eval_solve which would have basically all the solve algorithms. I agree that some kind of dispatching/pattern matching is useful, but I don't think the type of the expression is the correct thing to dispatch against.

@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented May 25, 2014

What are your thoughts on using a hint system, like dsolve? That would clearly separate the solve algorithms from their application.

@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented May 25, 2014

Yes, things like exp have an inverse method (I'm not clear why it's an instance method and not a class method, but there you have it).

# Maybe we can add the logic for lambert pattern here, better
# create a different function for it.
if g != S.One:
return [invt/g for invt in invert(h, symbol)]

This comment has been minimized.

@smichr

smichr May 26, 2014

Member

Although recursion is nice here, the current invert routine in solve uses a while loop which has its own elegance without the overhead of python function calls.

This comment has been minimized.

@smichr

smichr May 26, 2014

Member

My reference to invert is actually to _invert in solve.py.

# the conditionals greaterthan, lessthan are not defined
# for complex valued functions do something for it.
# also ploting for the piecewise functions doesn't work,
# it wll be easy to implement. Create a issue for it.

This comment has been minimized.

@smichr

smichr May 26, 2014

Member

If piecewise can be plotted then perhaps polygons can return a piece definition that plot_implicit could plot?

This comment has been minimized.

@asmeurer

asmeurer May 26, 2014

Member

I don't know why this comment is here.

Anyway, I don't know if plot_implicit supports Piecewise directly yet, but it does support Or, And, and conditionals, so in theory it should be doable by a logical translation using ITE.

This comment has been minimized.

@hargup

hargup May 26, 2014

Member

I don't know why this comment is here.

This is supposed to me note to myself. There will be more of them. I'll clean up such comments before the final rebase.

if g != S.One:
return [invt/g for invt in invert(h, symbol)]
if f.is_Add:
# f = g*h

This comment has been minimized.

@smichr

smichr May 26, 2014

Member

f = g + h

# f = g*h
g, h = f.as_independent(symbol)
if g != S.Zero:
return [invt - g for invt in invert(h, symbol)]

This comment has been minimized.

@smichr

smichr May 26, 2014

Member

return [invt.subs(symbol, symbol - g) for invt in invert(h, symbol)]?

This comment has been minimized.

@smichr

smichr May 26, 2014

Member

e.g. exp(x)+3 -> log(x - 3) not log(x) - 3

This comment has been minimized.

@smichr

smichr May 26, 2014

Member

The mul-block may also need work. solve_univariate(exp(x)*3-2,x) should give log((x + 2)/3).

This comment has been minimized.

@hargup

hargup May 27, 2014

Member

Thanks, I should write stronger test cases.

@hargup hargup referenced this pull request Dec 4, 2014

Open

Unstable submodule #8343

@hargup hargup changed the title [WIP] Rewriting the Univariate Solver Rewriting the Univariate Solver Dec 9, 2014

@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented Dec 9, 2014

So would you be OK with this if we removed the functions from __init__.py?

@hargup

This comment has been minimized.

Copy link
Member

hargup commented Dec 10, 2014

So would you be OK with this if we removed the functions from __init__.py?

I'm ok with it, @skirpichev ?

@skirpichev

This comment has been minimized.

Copy link
Contributor

skirpichev commented Dec 10, 2014

So would you be OK with this if we removed the functions from init.py?

Would we consider this as a public API?

@hargup

This comment has been minimized.

Copy link
Member

hargup commented Dec 11, 2014

@skirpichev

I feel if this PR doesn't get merged in near future then it will never be merged.
I'm involved in some other projects and I'm finding it difficult to contribute
enough time to sympy. Discussion and implementation of other new API might take
months and I'm unwilling to take part in it. And the chances that someone will
take up my work and modify it to get it merged are also slim. The
opportunity cost of not merging the PR at this stage is that the users will
potentially be forever stuck with the current messy solve. But we don't just
merge "anything". So I'll answer your question.

Would we consider this as a public API?

I don't know. But the question is, does it handle your concerns? Your major
concern is:

apparently we introduce public interface, which eventually will be broken.

I think you being too risk averse. Unlike incomplete code coverage or
objectively wrong results, the problem of an unstable API is not a very
concrete one. No API is ever perfect, there are always some cases which it
will not handle properly and there will always be cases to which it will not
generalize. Take the current solve for an example. It is messy at many
different levels but still there are use cases where it consistently gives the
needed result. Despite all its short comings people find it useful, and hence
people use it. If we wait for the perfectly stable API no new API will ever get
merged, even the current solve shouldn't have been merged at the first place.
And for the API which did get merged we weren't careful enough to see its
shortcomings. So the question is it is good enough for the users to find it
useful? And does it consistently solves a subset of the problem we are trying
to target? If the answer to both of the questions is yes then the apparent
unstable API shouldn't be a problem.

If you are still opposed to merging this PR then please suggest a clear and tangible
way forward.

@skirpichev

This comment has been minimized.

Copy link
Contributor

skirpichev commented Dec 11, 2014

On Thu, Dec 11, 2014 at 02:55:56AM -0800, Harsh Gupta wrote:

But the question is, does it handle your concerns?

If this still be a public API - No.

If you are still opposed to merging this PR then please suggest a clear
and tangible way forward.

I already did this.
#7523 (comment)

This is perfectly fine as a new experimental module.

@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented Dec 11, 2014

No, code that isn't exposed via __init__.py isn't public API. People are still free to use it, but it's "subject to change". The SymPy parsing code is another example of this. Ideally, we will refactor the parsers and it will all change, but currently quite a few people use it, because it's useful and the only way to do several things.

@skirpichev

This comment has been minimized.

Copy link
Contributor

skirpichev commented Dec 11, 2014

On Thu, Dec 11, 2014 at 02:44:10PM -0800, Aaron Meurer wrote:

No, code that isn't exposed via init.py isn't public API.

Well, then I'm ok with your suggestion.

@hargup hargup force-pushed the hargup:solver_refactor branch from 34e078a to 4f6e06a Dec 12, 2014

hargup added a commit to hargup/sympy that referenced this pull request Dec 12, 2014

Create the solveset module
The commit implements the solveset module in an attempt to
fix the problems with the current solve function. The commit
roughly follows my GSoC proposal
https://github.com/sympy/sympy/wiki/GSoC-2014-Application-Harsh-Gupta:-Solvers
and more details can be acquired through the discussion on the github pull
request sympy#7523

@hargup hargup force-pushed the hargup:solver_refactor branch from 4f6e06a to 0ec82e8 Dec 15, 2014

hargup added a commit to hargup/sympy that referenced this pull request Dec 15, 2014

Create the solveset module
The commit implements the solveset module in an attempt to
fix the problems with the current solve function. The commit
roughly follows my GSoC proposal
https://github.com/sympy/sympy/wiki/GSoC-2014-Application-Harsh-Gupta:-Solvers
and more details can be acquired through the discussion on the github pull
request sympy#7523

hargup added some commits Jun 22, 2014

Create the solveset module
The commit implements the solveset module in an attempt to
fix the problems with the current solve function. The commit
roughly follows my GSoC proposal
https://github.com/sympy/sympy/wiki/GSoC-2014-Application-Harsh-Gupta:-Solvers
and more details can be acquired through the discussion on the github pull
request #7523
@hargup

This comment has been minimized.

Copy link
Member

hargup commented Dec 15, 2014

Well, then I'm ok with your suggestion.

cool :) , I have removed solveset from the __init__.py file and also removed the @public decorators in the solveset file.

@skirpichev

This comment has been minimized.

Copy link
Contributor

skirpichev commented Dec 15, 2014

We have timeouts in test_solveset.py

@hargup hargup force-pushed the hargup:solver_refactor branch from 0ec82e8 to 641908e Dec 15, 2014

@hargup

This comment has been minimized.

Copy link
Member

hargup commented Dec 16, 2014

We have timeouts in test_solveset.py

fixed now, "The Travis CI build passed".

skirpichev added a commit that referenced this pull request Dec 16, 2014

Merge pull request #7523 from hargup/solver_refactor
Rewriting the Univariate Solver

@skirpichev skirpichev merged commit f28dd29 into sympy:master Dec 16, 2014

1 check passed

continuous-integration/travis-ci The Travis CI build passed
Details
@skirpichev

This comment has been minimized.

Copy link
Contributor

skirpichev commented Dec 16, 2014

Ok, this is in. Thank you, @hargup.

@hargup

This comment has been minimized.

Copy link
Member

hargup commented Dec 17, 2014

Thanks @mrocklin @skirpichev @asmeurer and @smichr for helping me out.

@mrocklin

This comment has been minimized.

Copy link
Member

mrocklin commented Dec 18, 2014

Hooray! Thank you Harsh for staying with this.
On Dec 16, 2014 11:06 PM, "Harsh Gupta" notifications@github.com wrote:

Thanks @mrocklin https://github.com/mrocklin @skirpichev
https://github.com/skirpichev @asmeurer https://github.com/asmeurer
and @smichr https://github.com/smichr for helping me out.


Reply to this email directly or view it on GitHub
#7523 (comment).

@skirpichev

This comment has been minimized.

Copy link
Contributor

skirpichev commented Dec 18, 2014

We still have a number of todo stuff. @hargup, could you look over this pr and fill bugreports? I'll try to help you asap.

@hargup

This comment has been minimized.

Copy link
Member

hargup commented Dec 18, 2014

We still have a number of todo stuff. @hargup, could you look over this pr and fill bugreports? I'll try to help you asap.

sure, We also have to add solveset to the docs, and a reference of solveset docs in the docs of solve.

@hargup hargup referenced this pull request Jan 1, 2015

Closed

Add docs for solveset #8725

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment