Skip to content
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

lambert solver for transolve #14972

Open
wants to merge 14 commits into
base: master
Choose a base branch
from

Conversation

Yathartha22
Copy link
Member

@Yathartha22 Yathartha22 commented Jul 25, 2018

This PR is built over #14792 and implements the lambert solver for _transolve in solveset. This has a dependency on bivariate.py for getting solutions for lambert type equations.

TODO's:

  • added lambert solver (_solve_lambert)

  • included tests that _solve_lambert solves, rest are added as XFAIL in test_solve_bivariate

  • include bivariate solver (bivariate_type)

  • Documentation and tests for solve_as_lambert

  • XFAIL (test_solve_bivariate) (will be done in a separate PR)

  • use solveset instead of solve in _lambert() in bivariate.py (will be done in a separate PR)

@aktech @smichr

Release Notes

  • solvers
    • added a new solver for lambert type equations as part of _transolve in solvers.solveset

@sympy-bot
Copy link

sympy-bot commented Jul 25, 2018

Hi, I am the SymPy bot (v158). I'm here to help you write a release notes entry. Please read the guide on how to write release notes.

Your release notes are in good order.

Here is what the release notes will look like:

  • solvers
    • added a new solver for lambert type equations as part of _transolve in solvers.solveset (#14972 by @Yathartha22)

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.6.

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it.

Click here to see the pull request description that was parsed.

This PR is built over #14792 and implements the lambert solver for `_transolve` in `solveset`. This has a dependency on `bivariate.py` for getting solutions for lambert type equations.

TODO's:

- [x] added `lambert` solver (`_solve_lambert`)

- [X] included tests that `_solve_lambert` solves, rest are added as XFAIL in `test_solve_bivariate`

- [x]  include `bivariate` solver (`bivariate_type`)

- [x] Documentation and tests for `solve_as_lambert` 

- XFAIL (`test_solve_bivariate`) (will be done in a separate PR)

- use `solveset` instead of `solve` in `_lambert()` in `bivariate.py` (will be done in a separate PR)

@aktech @smichr 

#### Release Notes

<!-- Write the release notes for this release below. See
https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
on how to write release notes. If there is no release notes entry for this PR,
write "NO ENTRY". The bot will check your release notes automatically to see
if they are formatted correctly. -->

<!-- BEGIN RELEASE NOTES -->
* solvers
  * added a new solver for lambert type equations as part of `_transolve` in `solvers.solveset`
<!-- END RELEASE NOTES -->

sympy/solvers/solveset.py Outdated Show resolved Hide resolved
sympy/solvers/solveset.py Outdated Show resolved Hide resolved
@Yathartha22
Copy link
Member Author

@smichr I have made the changes in #14792 once you approve there I will rebase it here.

@@ -1500,6 +1500,11 @@ def add_type(lhs, rhs, symbol, domain):

if new_eq is not None:
result = _solveset(new_eq, symbol, domain)
else:
# try solving as lambert before returning the result
ans = _solve_as_lambert(Eq(lhs, rhs), symbol)
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this the right place to invoke _solve_as_lambert? @smichr

@@ -1923,6 +1910,23 @@ def test_solve_lambert():
assert solveset_real(3*x + log(4*x), x) == \
FiniteSet(LambertW(Rational(3, 4))/3)


@XFAIL
def test_other_solve_lambert_():
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Added these tests as XFAIL because they are either solved as Mul or Pow

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I get no solution:

>>> solveset_real(-a*x + 2*x*log(x), x)
ConditionSet(x, Eq(-a*x + 2*x*log(x), 0), Reals)
>>> solveset_real(log(log(x - 3)) + log(x-3), x)
ConditionSet(x, Eq(x*log(x - 3) - 3*log(x - 3) - 1, 0), Reals)

@Yathartha22
Copy link
Member Author

Yathartha22 commented Aug 1, 2018

Some equations are solved by _solve_lambert after posifying the symbol, as in _tsolve. A few examples: x**2 - 2**x, (a/x + exp(x/2)).diff(x), x*log(x) + 3*x + 1, (a/x + exp(x/2)).diff(x, 2). Should we consider it in solveset too? What is the reason (maybe a proof) that making symbol positive would solve a large group of equations?


if lhs.is_Add:
result = add_type(equation, symbol, domain)
result = add_type(lhs, rhs, symbol, domain)
elif _is_lambert(lhs, symbol):
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

.is_Mul and .is_Pow here would become specific to lambert types, therefore a direct identification of lambert type is done here

@Yathartha22 Yathartha22 mentioned this pull request Aug 10, 2018
3 tasks
sympy/solvers/solveset.py Outdated Show resolved Hide resolved
@smichr
Copy link
Member

smichr commented Aug 10, 2018

What is the reason (maybe a proof) that making symbol positive would solve a large group of equations?

My guess is that this is so log transformations of powers would work.

sympy/solvers/solveset.py Outdated Show resolved Hide resolved
@Yathartha22
Copy link
Member Author

My guess is that this is so log transformations of powers would work.

Should solveset use this thing too or should we look for any other way?

@smichr
Copy link
Member

smichr commented Aug 12, 2018

Should solveset use this thing too or should we look for any other way?

the GSOC work suggests how solveset should handle this: make conditions part of what is returned. If all the conditions are met, the ConditionSet collapses, otherwise it is a sentinal for the requirements for a given solution.

If you don't use the log transformations I suspect you will have more cases to consider as you figure out what can and can't be solved by lambert.

…ustments in calling lambert and bivariate solver;
@Yathartha22
Copy link
Member Author

Not sure about this, firstly about raising NotImplementedError, it doesn't seem to go along with the traditional set output of solveset. Secondly as of now solveset is good enough in differentiating the type of equation, Lambert or bivariate, even has a way to tackle if not solved in the very first attempt (by making symbol positive), returns ConditionSet for the rest, all these things seems hard (as of now) to achieve with this way. Also _solve_lambert raises NotImplementedError, not sure how an exception is handled within an exception :-D.

Why not just add, a line in the documentation: mentioning about it returning partial solutions like in tsolve.

@jmig5776
Copy link
Member

I had tried to fix this. Please see this git diff @Yathartha22 @smichr

diff --git a/sympy/solvers/bivariate.py b/sympy/solvers/bivariate.py
index 51f4a9703f..3dca574e0c 100644
--- a/sympy/solvers/bivariate.py
+++ b/sympy/solvers/bivariate.py
@@ -142,6 +142,11 @@ def _lambert(eq, x):
         return []  # violated assumptions
     logarg = mainlog.args[0]
     b, c, X1 = _linab(logarg, x)
+    check_abs_root = False
+    if X1.is_Pow and X1.base is x and c == 0:
+        a = a*X1.exp
+        X1 = x
+        check_abs_root = True
     if X1 != X2:
         return []  # violated assumptions
 
@@ -152,14 +157,23 @@ def _lambert(eq, x):
         l = LambertW(d/(a*b)*exp(c*d/a/b)*exp(-f/a), k)
         # if W's arg is between -1/e and 0 there is
         # a -1 branch real solution, too.
-        if k and not l.is_real:
-            continue
+
+
         rhs = -c/b + (a/d)*l
 
         solns = solve(X1 - u, x)
         for i, tmp in enumerate(solns):
             solns[i] = tmp.subs(u, rhs)
-            sol.append(solns[i])
+            if solns[i].is_real:
+                sol.append(solns[i])
+
+        if check_abs_root:
+            l = LambertW(-d/(a*b)*exp(c*d/a/b)*exp(-f/a), k)
+            rhs = -c/b + (a/d)*l
+            solns2 = solve(X1 - u, x)
+                solns2[i] = tmp.subs(u, rhs)
+                if solns2[i].is_real:
+                    sol.append(solns2[i])
     return sol
 
 
@@ -206,8 +220,8 @@ def _solve_lambert(f, symbol, gens):
         raise NotImplementedError()
 
     if lhs.is_Mul:
-        lhs = expand_log(log(lhs))
-        rhs = log(rhs)
+        lhs = expand_log(log(abs(lhs)))
+        rhs = log(abs(rhs))
 
     lhs = factor(lhs, deep=True, fraction=False)
     # make sure we have inverted as completely as possible
@@ -235,16 +249,16 @@ def _solve_lambert(f, symbol, gens):
         mainlog = _mostfunc(lhs, log, symbol)
         if mainlog:
             if lhs.is_Mul and rhs != 0:
-                soln = _lambert(log(lhs) - log(rhs), symbol)
+                soln = _lambert(log(abs(lhs)) - log(abs(rhs)), symbol)
             elif lhs.is_Add:
                 other = lhs.subs(mainlog, 0)
                 if other and not other.is_Add and [
                         tmp for tmp in other.atoms(Pow)
                         if symbol in tmp.free_symbols]:
                     if not rhs:
-                        diff = log(other) - log(other - lhs)
+                        diff = log(abs(other)) - log(abs(other - lhs))
                     else:
-                        diff = log(lhs - other) - log(rhs - other)
+                        diff = log(abs(lhs - other)) - log(abs(rhs - other))
                     soln = _lambert(expand_log(diff), symbol)
                 else:
                     #it's ready to go
@@ -267,7 +281,7 @@ def _solve_lambert(f, symbol, gens):
         if mainexp:
             lhs = collect(lhs, mainexp)
             if lhs.is_Mul and rhs != 0:
-                soln = _lambert(expand_log(log(lhs) - log(rhs)), symbol)
+                soln = _lambert(expand_log(log(abs(lhs)) - log(abs(rhs))), symbol)
             elif lhs.is_Add:
                 # move all but mainexp-containing term to rhs
                 other = lhs.subs(mainexp, 0)
@@ -277,7 +291,7 @@ def _solve_lambert(f, symbol, gens):
                     rhs.could_extract_minus_sign()):
                     mainterm *= -1
                     rhs *= -1
-                diff = log(mainterm) - log(rhs)
+                diff = log(abs(mainterm)) - log(abs(rhs))
                 soln = _lambert(expand_log(diff), symbol)
 
     # 3) d*p**(a*B + b) + c*B = R
@@ -289,13 +303,13 @@ def _solve_lambert(f, symbol, gens):
         if mainpow and symbol in mainpow.exp.free_symbols:
             lhs = collect(lhs, mainpow)
             if lhs.is_Mul and rhs != 0:
-                soln = _lambert(expand_log(log(lhs) - log(rhs)), symbol)
+                soln = _lambert(expand_log(log(abs(lhs)) - log(abs(rhs))), symbol)
             elif lhs.is_Add:
                 # move all but mainpow-containing term to rhs
                 other = lhs.subs(mainpow, 0)
                 mainterm = lhs - other
                 rhs = rhs - other
-                diff = log(mainterm) - log(rhs)
+                diff = log(abs(mainterm)) - log(abs(rhs))
                 soln = _lambert(expand_log(diff), symbol)
 
     if not soln:

This gives the required 3 solutions in (solveset_real((1/x + exp(x/2)).diff(x), x))

>>> (solveset_real((1/x + exp(x/2)).diff(x), x))
{4*LambertW(-sqrt(2)/4), 4*LambertW(sqrt(2)/4), 4*LambertW(-sqrt(2)/4, -1)}

>>> N((solveset_real((1/x + exp(x/2)).diff(x), x)))
{−5.23573322613363,−2.97592413099635,1.07967055380564}

Although I am trying to fix it generally. Some tests are still failing due to introduction of abs. It might take some time. Please give me suggestions about above approach.

@jmig5776
Copy link
Member

jmig5776 commented Mar 23, 2019

Basically the problem here we are tackling is that _solve_lambert tries to convert all type of lambert expressions to this form F(X, a..f) = a*log(b*X + c) + d*X + f = 0 in which while converting them to log some solutions are missed like the following case:

solveset_real((1/x + exp(x/2)).diff(x), x)

The equation that goes in _lambert here is like this

_r + 2*log(_r**2) - 2*log(2)

So our job is to get the 2 inside of log to outside to convert to the required form otherwise it is returning ConditionSet.
In equation like this

solveset_real((a/x + exp(x/2)).diff(x, 2), x)

The equation that goes in _lambert in my above git diff is

_r + 2*log(_r**2*Abs(_r)) - 2*log(8*Abs(a))

Here I dont know how to get 3(2+1 of _r and abs(_r)) outside to get into this form. Please help here @smichr.

And the point regarding domain. There will be no need to change the symbol to positive if we are able to convert the equation in real form to required form in lambert (in a proper way of considering abs). Here changing symbol to positive is passing most of the test cases because they are designed in such way.

This is I think the problem and above is my plan. Please review @Yathartha22 , @smichr

@smichr
Copy link
Member

smichr commented Mar 24, 2019

Why not just add, a line in the documentation

If it is reported by solveset it has to adhere to the stated claim:

    Set
        A set of values for `symbol` for which `f` is True or is equal to
        zero. An `EmptySet` is returned if `f` is False or nonzero.
        A `ConditionSet` is returned as unsolved object if algorithms
        to evaluate complete solution are not yet implemented.

    `solveset` claims to be complete in the solution set that it returns.

That's why I feel this is not ready to add to solveset. Can you return the union of what is solved and a ConditionSet containing the equation?

What is hindering being able to return all solutions?

@smichr
Copy link
Member

smichr commented Mar 24, 2019

to get into this form

>>> logcombine(log(x)+log(2)+3)
log(2*x) + 3

logcombine won't combine arguments if the sign is going to be unknown.

from sympy import *
from sympy.abc import x,y
print(logcombine(y + 2*log(y**2*Abs(y)) - 2*log(8*Abs(x))))
print(logcombine(log(x)+log(2)))
print(logcombine(log(x**2)+log(2)))
print(logcombine(log(x**3)+log(2)))
var('x y',positive=1)
print(logcombine(y + 2*log(y**2*Abs(y)) - 2*log(8*Abs(x))))
print(logcombine(log(x)+log(2)))
print(logcombine(log(x**2)+log(2)))
print(logcombine(log(x**3)+log(2)))
y + 2*log(y**2*Abs(y)) - 2*log(8*Abs(x))
log(x) + log(2)
log(x**2) + log(2)
log(x**3) + log(2)
y + log(y**6/(64*x**2))
log(2*x)
log(2*x**2)
log(2*x**3)

It should combine a positive with an unknown, however, since this would not change the sign. See #6770.

BTW, the condition that an argument be positive is something that could get back with the solution as part of the condition set. If the conditions are met, they won't be retained; if they aren't, then the user will know the requirements for the solution.

@jmig5776
Copy link
Member

jmig5776 commented Apr 7, 2019

@Yathartha22 May I cherry pick your the commit of your this PR and open a new PR to complete lambert? Let me know ASAP.

@Yathartha22
Copy link
Member Author

Can you return the union of what is solved and a ConditionSet containing the equation?

Yeah I guess. I will have a look at it.

May I cherry pick your the commit ...

I don't think it would be a good idea, because the issue with the PR is not returning all solutions, which possibly can be achieved by improving solve_lambert as you were mentioning above. So it would be better you open a separate PR for only this purpose (why to unnecessary increase the diff by cherry picking).
Here we can return union of solutions with conditionset and once the solve_lambert is improved solveset will automatically include all those solutions.

@jmig5776
Copy link
Member

jmig5776 commented Apr 8, 2019

Okk @Yathartha22 , I will be opening a PR to improve solve_lambert soon.

@czgdp1807
Copy link
Member

@Yathartha22 I believe that @jmig5776 has made PR to continue this one and it has been merged. So can we close this?

@Yathartha22
Copy link
Member Author

The changes done by Jogi has been in the bivariate which will help to get the missing solutions for this PR and then it can be merged. I think final touchup need to be done in this PR.
I have been busy these days, I will try to finish whatever is left with this PR soon.

@MaqOwais
Copy link
Contributor

I've gone through all the work done by @jmig5776, @smichr and other's related to lambert solver where , further we have to extend it to _transolve.

I'm really glad to work on it ,...Please support me if I go wrong.

@czgdp1807
Copy link
Member

Closing in favour of #18759

@czgdp1807 czgdp1807 closed this Mar 12, 2020
@Yathartha22
Copy link
Member Author

@czgdp1807 I don't think this PR should be closed in favour of #18759 as the latter doesn't seem to be convincing also this PR has a lot of work and important commits.
And why is it that I am unable to reopen my own PR, that's weird.

@oscarbenjamin oscarbenjamin reopened this Apr 18, 2020
@oscarbenjamin
Copy link
Collaborator

@Yathartha22 github shows you as a "contributor" rather than a "member". I think that a non-member can not reopen a PR/issue that was closed by a member.

If you were previously a member then it's possible that it has automatically expired. I think this happens after 1 year if you don't commit anything but you can email @asmeurer and he will make you a member again.

@Yathartha22
Copy link
Member Author

Thanks @oscarbenjamin for the clarification and reopening the PR, I will have a look at it.

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

Successfully merging this pull request may close these issues.

8 participants