log solver for transolve

[Yathartha22]

This PR is an extension to #14736 and is built over its branch.

The PR aims to implement log solver as part of transolve.

TODO's:

• Documentation for _solve_ and _is_ routines for log

• Tests for the helpers

• Improvement in _solve_ routine to handle corner cases.

Release Notes

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

[Yathartha22]

I am not sure if this is the correct way to get some features from a different PR that this branch is dependent on.

ping @aktech @smichr @asmeurer

[smichr]

I am not sure if this is the correct way to get some features from a different PR that this branch is dependent on.

If the other branch is very stable, this is ok...but I have often run into rebase issues when doing this.

[Yathartha22]

If the other branch is very stable, this is ok...but I have often run into rebase issues when doing this.

I guess it won't be stable, it will have few commits more. I will give it a try though otherwise I will wait for main branch to get merge and add a new PR.

[smichr]

I believe you do not want to do this. You want to do the solution and domain checking all in one place, not duplicate it in every helper for every type of equation.

[Yathartha22]

Will be adding documentation in the next commit.

Logic has been changed to check each term of the expression of log form rather than previously used logcombine method:

• using logcombine can be risky as it may manipulate the equation like log(x) - log(2*x) would reduce to log(1/2) and the routine would return False as there will be no symbol.

[Yathartha22]

This is for expressions like log(x)**2.

[Yathartha22]

This check is done so as to remove unwanted solutions form the original result.
This method is implemented rather than using checksol because there are some unwanted solution that would still creep into the actual result.

Take for example log(3*x) - log(-x + 1) - log(4*x + 1), the result in general is -1/2 and 1/2 but the negative solution is only valid if we consier complex valued logairthm but since we are concerned with the Real domain as of now only natural logairthm will work which is undefined for negative as well as for 0 therefore -1/2 can't be the solution in real doamin.
If the result is verified using checksol the subsitution would cause the cancelling of I*pi and hence it would return True, making -1/2 as the solution

In this method each solution is tried in each term of the expression, if its result is not in a real form then it is not the solution (This is specifically beneficial for checking for logarithmic expressions).

[Yathartha22]

I will add the remaining tests as well as the documentation once we agree on the current changes.

ping @aktech @smichr

[Yathartha22]

Will be adding docstring in the next commit.

[Yathartha22]

@smichr @aktech I have started a discussion to arrive at a conclusion for logarithmic conclusion. Have a look.

[Yathartha22]

• I have removed checking log singularities as of now unless we get to a conclusion here. This needs more discussion.

• I have included the use of flag variable (and this time it is explicitly declared). Some equations that do not get solved by exponential or logarithmic solver but are still one of these (eg: equations solvable as lambert) can get into a recursion.

ping @aktech @smichr

[smichr]

->containing (sp)

[smichr]

suggest combining into for expr_arg in _term_factors(f):

[Yathartha22]

Just used to improve readability. One could at once get what is returned from _term_factors. Should I change it?

[smichr]

is the helper necessary? if isinstance(expr_arg, log) and symbol in expr_arg.free_symbols: will be less than 72 characters.

[smichr]

oh...you use it again below. Then how about prefacing with

base = expr_arg.base if isinstance(expr_arg, Pow) else expr_arg
if isinstance(base, log) and symbol in base.free_symbols:
    return True

[Yathartha22]

This can be done. Thanks!

[smichr]

Like exponential equations, you should describe the definition of being logarithmic.

[Yathartha22]

Seems that calling _solveset(...) within Piecewise is not good enough, because irrespective of the condition the function is invoked. (I have asked about this scenario in gitter)
I guess using if won't make a difference.

ping @aktech @smichr

[smichr]

My main concern is that if there are conditions that are neither false nor true -- just undecided -- I would like to see the solution and the conditions come back, not a generic Eq(foo, 0). Then a person doesn't have to go through the solving process when parameters are ready to be substituted. Give the maximum information for their query. (I'll see what I can make of the Piecewise problem.)

[smichr]

This will handle the recursion error:

old = [log(a), log(-b)]
new = list(expand_log, old))
if new == old:
    return unsolved_result
return _solveset(new[0] - new[1], symbol, domain)

Nothing in solveset is ready to handle Piecwise solutions. It expects sets to be returned. That being the case, I don't know how the solution with conditions can be returned. You might as well just return functions (instead of solutions) from the helpers. In this case, new[0] - new[1].

One possible use of conditions, however, would be to return EmptySet if the conditions are False.

[Yathartha22]

That being the case, I don't know how the solution with conditions can be returned.

IMO we can do this in two ways:

• Raise a ValueError, with a message like: "Proper assumptions needs to be givien to solve the given type of exponential equation "
• Force the equation to get solved but return a ConditionSet (a different than ususal) with the required assumption, like:
  ConditionSet(x, And(a_base>0, b_base>0, Eq(im(a_exp), 0), Eq((im(b_exp), 0))), {solutions})

for equation a**x - b**x (vanilla), we can have
ConditionSet(x, And(a>0, b>0, Eq(im(x), 0), Eq((im(x), 0))), {0})

One possible use of conditions, however, would be to return EmptySet if the conditions are False.

Returning EmptySet could be misleading in a way as it would be pointing out to two things: either x doesn't have a solution in the respective domain or the assumptions are invalid.

[smichr]

doesn't have a solution in the respective domain or the assumptions are invalid.

In either case it seems the appropriate answer is EmptySet. e.g. that there is no answer in a domain doesn't mean there isn't an answer. If the user gives a domain and we know the answer is not in that domain then the empty set is appropriate:

>>> solveset(x**3+8,x,Interval(0,oo))
EmptySet()

Regarding the condition set return value, I think you are right. What about:

return `Intersection(ConditionSet(...), domain)`

[Yathartha22]

return Intersection(ConditionSet(...), domain)

Will Intersection make any difference? because I think solutions will be intersected with the domain by _solveset itself

[Yathartha22]

Just one query about this scenario. Forcing the equation to be solved prompts solveset to return equation as an intersection with the domain.
For eg: if I force the following equation to get solved

>>> a, b, x, y = symbols('a b x y')
>>> solveset(a**x - b**(x + 1), x, S.Reals)
Intersection(Reals, {log(b)/(log(a) - log(b))})

This is something similar to what we get in logarithmic eqautions. So is returning this kind of ConditionSet still appropriate or should we try force solving the equation (to get something as above)?

[smichr]

I think solutions will be intersected with the domain by _solveset itself

If this is the case then passing domain to the helpers is not necessary. Once can helper-solve in the complex domain and _solveset will apply the limiting domain in the intersection.(?)

[Yathartha22]

Once can helper-solve in the complex domain and _solveset will apply the limiting domain in the intersection.(?)

But still, we will be needing to pass the "domain" parameter in _solveset (S.Complexes is not the default argument for domain).

[Yathartha22]

@smichr @aktech I am hoping to get the PR merge before this week ends (last week for GSoC). It would great that we make it ready by then. I can open another PR for any other related issues (that is not breaking this PR)

[aktech]

@Yathartha22 I am fine with merging this as long as @smichr doesn't have any concerns.

[smichr]

diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index 83b6ffb..1626d8d 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -1113,10 +1113,9 @@ def _solve_exponential(lhs, rhs, symbol, domain):
         Eq(im(b_exp), 0)
         )
 
-    if conditions is S.true:
-        return _solveset(expand_log(log(a)) - expand_log(log(-b)), symbol, domain)
-
-    return unsolved_result
+    neweq = expand_log(log(a), force=True) - expand_log(log(-b), force=True)
+    soln = _solveset(neweq, symbol, domain)
+    return ConditionSet(symbol, conditions, soln)
 
 
 def _is_exponential(f, symbol):
@@ -1528,9 +1527,6 @@ def add_type(lhs, rhs, symbol, domain):
     else:
         result = rhs_s
 
-    if isinstance(result, ConditionSet):
-        result = unsolved_result
-
     return result
 
 
diff --git a/sympy/solvers/tests/test_solveset.py b/sympy/solvers/tests/test_solveset.py
index 37c9166..a6879cf 100644
--- a/sympy/solvers/tests/test_solveset.py
+++ b/sympy/solvers/tests/test_solveset.py
@@ -1867,8 +1867,10 @@ def test_exponential_symbols():
         (x - z)/log(2)) - FiniteSet(0)
 
     a, b, x, y = symbols('a b x y')
+    assert solveset(a**x - b**x, x) == ConditionSet(
+        x, (a > 0) & (b > 0) & Eq(im(x), 0), {0})
     assert solveset_real(a**x - b**x, x) == ConditionSet(
-        x, Eq(a**x - b**x, 0), S.Reals)
+        x, (a > 0) & (b > 0), {0})
 
 
 @XFAIL

The docstrings needs a touch up if you use this diff to reflect that ConditionSet may be returned if there are conditions that must be met for the solution to be valid. It seems a little ill-defined right now as to who decides there is no solution, and whether the difference between 'can't solve' and 'no solution' is clearly being made.

If used, a few tests need editing, too, to put the constant on the rhs of the Eq:

___________ sympy\solvers\tests\test_solveset.py:test_conditionset ____________
Traceback (most recent call last):
  File "c:\users\leslie\sympy\sympy\solvers\tests\test_solveset.py", line 954, i
n test_conditionset
    ) == ConditionSet(x, Eq(x**2 + x*sin(x) - 1, 0), S.Reals)
AssertionError
_______________________________________________________________________________
_________ sympy\solvers\tests\test_solveset.py:test_expo_conditionset _________
Traceback (most recent call last):
  File "c:\users\leslie\sympy\sympy\solvers\tests\test_solveset.py", line 1843,
in test_expo_conditionset
    x, Eq(2**x - exp(x) - 3, 0), S.Reals)
AssertionError

[Yathartha22]

Maybe we define the "conditions " for real domain only:

conditions=True
if domain.is_subset(Reals):
        conditions = And(
            a_base > 0,
            b_base > 0,
            Eq(im(a_exp), 0),
            Eq(im(b_exp), 0))

and solveset(a**x - b**x, x, S.Complexes) should return simply {0}.

Note: Even though doing this we will be getting the solutions in Complex domain as well (but they won't be general ones). I would rather make a separate PR for solving exponential equations in complex domain (genral solutions)

who decides there is no solution, and whether the difference between 'can't solve' and 'no solution' is clearly be made.

I guess _solveset takes care of this

[Yathartha22]

if used, a few tests need editing, too, to put the constant on the rhs of the Eq:

@smichr I have already made the changes in the latest commit below this comment. I guess the branch you are using is not updated.

[Yathartha22]

This will be the behaviour

>>> f = a**x - b**x
>>> solveset(f, x, S.Reals)
ConditionSet(x, (a > 0) & (b > 0), {0})
>>> solveset(f, x, S.Complexes)
{0}
>>> a, b = symbols('a b', positive=True)
>>> solveset(f, x, S.Reals)
{0}

[smichr]

This will be the behaviour

There are still assumptions when the domain is complex:

>>> q=a**x+b**(x+1)
>>> solveset(q,x)
{(log(b) + I*pi)/(log(a) - log(b))}

We should require that a and b be nonzero and that a not equal b: Ne(a, b), Ne(a, 0), Ne(b, 0)

[Yathartha22]

Oh yes I missed this for complex domain. Thanks!
I guess Ne(a, b) can be used for real domain as well?
I see a problem with Ne(a, b) for equations having exp as their base: -9*exp(-2*x + 5) + 4*exp(3*x + 1). Should we leave it?

[smichr]

Please sharpen a pencil and check out the assertions in the following code comments:

diff --git a/sympy/solvers/solveset.py b/sympy/solvers/solveset.py
index 5ac6781..15f4849 100644
--- a/sympy/solvers/solveset.py
+++ b/sympy/solvers/solveset.py
@@ -13,7 +13,7 @@
 from __future__ import print_function, division
 
 from sympy.core.sympify import sympify
-from sympy.core import (S, Pow, Dummy, pi, Expr, Wild, Mul, Equality,
+from sympy.core import (S, Pow, Dummy, pi, Expr, Wild, Mul,
                         Add)
 from sympy.core.containers import Tuple
 from sympy.core.facts import InconsistentAssumptions
@@ -27,10 +27,11 @@
 from sympy.functions import (log, Abs, tan, cot, sin, cos, sec, csc, exp,
                              acos, asin, acsc, asec, arg,
                              piecewise_fold, Piecewise)
+from sympy.functions.elementary.complexes import sign, im
 from sympy.functions.elementary.trigonometric import (TrigonometricFunction,
                                                       HyperbolicFunction)
 from sympy.functions.elementary.miscellaneous import real_root
-from sympy.logic.boolalg import And
+from sympy.logic.boolalg import And, ITE
 from sympy.sets import (FiniteSet, EmptySet, imageset, Interval, Intersection,
                         Union, ConditionSet, ImageSet, Complement, Contains)
 from sympy.sets.sets import Set
@@ -246,10 +247,14 @@ def _invert_real(f, g_ys, symbol):
                 s, b = integer_log(rhs, base)
                 if b:
                     return _invert_real(expo, FiniteSet(s), symbol)
-            elif rhs is S.One:
-                # special case: 0**x - 1
-                return (expo, FiniteSet(0))
-            return (expo, S.EmptySet)
+            elif base.is_zero:
+                one = Eq(rhs, 1)
+                if one == True:
+                    # special case: 0**x - 1
+                    return (expo, FiniteSet(0))
+                elif one == False:
+                    return (expo, S.EmptySet)
+            return (f, g_ys)
 
 
     if isinstance(f, TrigonometricFunction):
@@ -877,7 +882,7 @@ def _solveset(f, symbol, domain, _check=False):
     inverter = lambda f, rhs, symbol: _invert(f, rhs, symbol, domain)
 
     result = EmptySet()
-
+    print(f)
     if f.expand().is_zero:
         return domain
     elif not f.has(symbol):
@@ -1104,27 +1109,64 @@ def _solve_exponential(lhs, rhs, symbol, domain):
         return unsolved_result
 
     a, b = list(ordered(lhs.args))
-    a_term = a.as_independent(symbol)[1]
-    b_term = b.as_independent(symbol)[1]
+    a_co, a_term = a.as_independent(symbol)
+    b_co, b_term = b.as_independent(symbol)
+
+    # if we don't know if the equation really is two args then
+    # we can't return a solution as though it were and there
+    # is no way to represent a Piecewise solution in solveset
+    # so we have to return an unsolved result
+    if And(Ne(a_co, 0), Ne(b_co, 0)) != True:
+        # indeterminate:
+        # if a==0 or c==0, the solution of a*b**f(x) + c*d**g(x)
+        # is empty set if d or b (respectively) is not zero else
+        # if a==0 and d==0 then g(x)=1 or if c==0 and b==0 then
+        # f(x) == 1; if a and c are 0 the solution is the domain
+        return unsolved_result
 
     a_base, a_exp = a_term.base, a_term.exp
     b_base, b_exp = b_term.base, b_term.exp
 
-    from sympy.functions.elementary.complexes import im
+    if ITE(
+            Eq(a_base, b_base),
+            Ne(a_exp, b_exp),
+            ITE(
+                Eq(a_exp, b_exp),
+                Ne(a_base, b_base),
+                Ne(a, b))) != True:
+        # indeterminate:
+        # if bases are equal there are two possible solutions:
+        #  if(b=d=B) then a*b**f(x) + c*d**g(x)
+        #  is a*B**f(x) + c*B**g(x); if the exponents are the same, f=g=h, then
+        #  the solution is given by B**h(x)*(a + c); if B is zero then the solution
+        #  is h(x) != 0 else the solution is the empty set.
+        # If the exponents are the same, there are 3 possible solutions
+        #  if f=g=h, then we have a*b**h(x)+c*d**h(x). The case of b==d
+        #  was handled above and there are two cases: h(x) != 0 or the empty set.
+        # If b!=d then log(a) + h(x)*log(b) - log(-c) - h(x)*log(d)
+        # which gives log(a/-c) + h(x)*log(b/d) and h(x) = -log(a/-c)/log(b/d)
+        return unsolved_result
 
+    L, R = map(lambda i: expand_log(log(i), force=True), (a, -b))
     if domain.is_subset(S.Reals):
         conditions = And(
+            Ne(sign(a_co), sign(b_co)),
             a_base > 0,
             b_base > 0,
             Eq(im(a_exp), 0),
-            Eq(im(b_exp), 0))
+            Eq(im(b_exp), 0),
+            )
+        solutions = _solveset(L - R, symbol,
+            # when domain is Reals, conditions will take care
+            # of keeping the solution Real, otherwise we
+            # need the clipped solution from solveset
+            S.Complexes if domain is S.Reals else domain)
     else:
         conditions = And(
             Ne(a_base, 0),
-            Ne(b_base, 0))
-
-    log_type_equation = expand_log(log(a), force=True) - expand_log(log(-b), force=True)
-    solutions = _solveset(log_type_equation, symbol, domain)
+            Ne(b_base, 0),
+            )
+        solutions = _solveset(L - R, symbol, domain)
 
     return ConditionSet(symbol, conditions, solutions)
 
@@ -1884,7 +1926,7 @@ def linear_eq_to_matrix(equations, *symbols):
 
     for equation in equations:
         f = sympify(equation)
-        if isinstance(f, Equality):
+        if isinstance(f, Eq):
             f = f.lhs - f.rhs
 
         # Extract coeff of symbols

[smichr]

The comments are not "production ready" but I have to leave this to you while I end my own week's break and finish things delayed.

[smichr]

I opened an issue for the mishandling of a single exponential term, so that doesn't have to be fixed. I think this routine is giving an answer in cases where there is not enough information to return a Set solution, but that can be fixed later. The main thing is that the structure is in place. I think the GSOC has a category for "future improvements". So there's no reason to delay this PR.

Congrats. I will be checking during breaks today to see if there is anything else you need reviewed.

[Yathartha22]

Well, thanks @smichr for the merge and all the reviews :).

I was actually working on diff that you suggested.
Few points on it:

1. The generalization of the equations was great but what I deduced (and correct me if I am wrong) that these types are not necessaary to be handled here, because they can be handled in solveset (specially coefficient being 0 and the base being 0 cases).
2. condition like Ne(sign(a_co), sign(b_co) can be included in conditions.
3. The use of map (a better pythonic way)

For points 2, 3 and the issue you created, I can create a seperate PR addressing all these points.
Anything else that needs to be fixed from this branch?

I will be checking during breaks today to see if there is anything else you need reviewed.

Next aim would be adding a lambert solver #14972. I will rebase it with the current master.