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

transcendental equation solver for solveset: transolve #14736

Merged
merged 34 commits into from Jul 17, 2018

Conversation

Projects
None yet
5 participants
@Yathartha22
Contributor

Yathartha22 commented May 23, 2018

This PR aims at extending capability of solveset by implementing transolve that is a helper to solveset for solving transcendental equations.

The following things are done:

  • API design and flow of transolve
  • Exponential solver
  • import test from solve
  • Fix pep8
  • address comments

The following needs to be done:

  • solving exponential equation in complex domain~ (will be solved later, added tests as XFAIL as of now)
  • XFails

  • Documentation with proof of correctness.

this is an alternative to #13045.
#10864

ping @aktech @smichr

@Abdullahjavednesar Abdullahjavednesar requested review from aktech and smichr May 23, 2018

f2 = x**x
assert solveset(e1, x, S.Reals) == FiniteSet(-3*log(2)/(-2*log(3) + log(2)))
assert simplify(solveset(e2, x, S.Reals)) == FiniteSet(4/S(15))

This comment has been minimized.

@Yathartha22

Yathartha22 May 23, 2018

Contributor

Need to return a simplified result for such equations.
eg: 2**x - 32 currently gives {log(32)/log(2)}, should rather be {5}.
Is it a good idea to use simplify() here?

This comment has been minimized.

@asmeurer

asmeurer May 23, 2018

Member

Raw simplify can be a very expensive function call and is generally avoided in library code. In the old solve there is a flag to enable or disable simplify.

Is it possible instead to invert powers using a two argument log, like log(32, 2)?

This comment has been minimized.

@Yathartha22

Yathartha22 May 24, 2018

Contributor

@asmeurer it seems that two argument log does not work well with imageset. I tried

>>> g = FiniteSet(32)
>>> imageset(Lambda(n, log(n, 2)), g)
{log(32)/log(2)}
assert solveset(z**x - y, x, S.Reals) == Intersection(S.Reals, FiniteSet(log(y)/log(z)))
assert solveset(y - a*x**b, x) == FiniteSet((y/a)**(1/b))
assert solveset(Poly(exp(x) + exp(-x) - 4)) == FiniteSet(log(-sqrt(3) + 2), log(sqrt(3) + 2))

This comment has been minimized.

@Yathartha22

Yathartha22 May 23, 2018

Contributor

This equation can be handled by solve_decomposition.

This comment has been minimized.

@asmeurer

asmeurer May 23, 2018

Member

Why is Poly used here?

This comment has been minimized.

@Yathartha22

Yathartha22 May 23, 2018

Contributor

I think Poly should be removed as there is no condition in solveset that would extract expression from Poly as it does in solve. I might have missed it.
Can we have a check of such case just to increase the coverage?

if isinstance(f, Poly):
            f = f.as_expr()
rhs = -b
lhs = expand_log(log(lhs), force=True)
rhs = expand_log(log(rhs), force=True)

This comment has been minimized.

@asmeurer

asmeurer May 23, 2018

Member

force=True does operations that are not mathematically correct for all values. Is it possible for this function to return wrong results because of this?

This comment has been minimized.

@smichr

smichr May 26, 2018

Member

Yes, log(x*y) will expand into log(x)+log(y) but that is not always true. e.g. if x == y == -1, the former is 0 while the latter is 2*pi*I in SymPy.

@aktech

This comment has been minimized.

Member

aktech commented May 23, 2018

@Yathartha22 This PR doesn't fulfill the minimum requirements for review, which we discussed in the first meeting. I don't see this comment being taken care of.

rhs = rhs_s.args[0]
# do we need a loop for rhs_s?
# Can't determine a case as of now.
if lhs == symbol:

This comment has been minimized.

@Yathartha22

Yathartha22 May 24, 2018

Contributor

This case is added so that even if one tries to solve the equation by calling transolve, result can be obtained.
cases like: 2**x -32 will be handled here.

This comment has been minimized.

@smichr
@Yathartha22

This comment has been minimized.

Contributor

Yathartha22 commented May 25, 2018

ping @aktech

transolve's API and its flow is simple to understand. Unlike _tsolve which
computes by solving with lots of recursive calls.
transolve eases the task by reducing it to a two step procedure as
dicussed below.

This comment has been minimized.

@smichr

smichr May 26, 2018

Member

typos

inver_real - invert_real
eqaution - equation
generalised - generalized (I know, variant) :-)
dicussed - discussed
doamin - domain
eqquation - equation

{-log(375)/(-log(5) + 2*log(3))}
"""
if 'tsolve_saw' not in flags:

This comment has been minimized.

@smichr

smichr May 26, 2018

Member

flags.setdefault('tsolve_saw', [])

This comment has been minimized.

@smichr

smichr May 26, 2018

Member

Or all in one:

if f in flags.setdefault('tsolve_saw', []):
    return...
flags['tsolve_saw'].append(f)
else:
flags['tsolve_saw'].append(f)
inverter = invert_real if domain.is_subset(S.Reals) else invert_complex

This comment has been minimized.

@smichr

smichr May 26, 2018

Member

A call to invert_complex(..., domain) will handle which inverter to call (like calling to solveset will just work when you pass a domain).

This comment has been minimized.

@Yathartha22

Yathartha22 May 26, 2018

Contributor

Yeah, we can do that. Will add a comment so as to avoid confusion.

result = rhs_s
elif lhs.is_Add:
equation = factor(powdenest(lhs-rhs))

This comment has been minimized.

@smichr

smichr May 26, 2018

Member

you have to capture denominators containing possible singularities before doing lhs-rhs or else terms with denominators may cancel, e.g. x*y - 1/x = x*z - 1/x for which x = 0 is not a solution.

result = S.EmptySet
else:
rhs = rhs_s.args[0]
# do we need a loop for rhs_s?

This comment has been minimized.

@smichr

smichr May 26, 2018

Member
assert (rhs_s.args) == 1  # can't think of a case where there would be more than one

This comment has been minimized.

@aktech

aktech Jun 12, 2018

Member

+1

@Yathartha22 This comment doesn't seem to be addressed as of yet. Never assume anything, always handle all the cases.

This comment has been minimized.

@Yathartha22

Yathartha22 Jun 13, 2018

Contributor

I will add a loop here rather an assert statement.

result = rhs_s
elif lhs.is_Add:
equation = factor(powdenest(lhs-rhs))

This comment has been minimized.

@smichr

smichr May 26, 2018

Member

Would factor_terms be sufficient for what you are expecting? See its options, too, for handling radicals.

This comment has been minimized.

@Yathartha22

Yathartha22 May 26, 2018

Contributor

I will have a look over it. But here I am trying to create a P.O.S of the equation.

ii) Invoking the respective helper to solve the equation.
Once identified what family the equation belongs, respective

This comment has been minimized.

@smichr

smichr May 26, 2018

Member

Once the family type of the equation has been identified, the corresponding helper is called...

Helper function for solving exponential equations.
This function solves exponential equation by using logarithms.
Exponents are converted to a more general form of logs which can further be

This comment has been minimized.

@smichr

smichr May 26, 2018

Member

...which can be solved by solveset.

The design of `\_transolve` makes adding a new class of equation
solver easy.
Decide from where the *identification* *helper* will be invoked.

This comment has been minimized.

@smichr

smichr Jul 4, 2018

Member

Is there any transcendental equation (other than something like cos(x) = 0) that can be solved/handled by transolve that is not an Add?

This paragraph is confusing with the example put in the middle of the instructions. "Once done" and "add a call" are not clearly defined. It would be better to give a code snippet like

if eq.is_Add:
    ...
    if _is_expo(eq, x):
        new_eq = _solve_expo(eq, x)
    ...

Also, what to name the result should probably be specified. Are different names used for the case when a solution is returned vs the case where a rewritten form is returned?

This comment has been minimized.

@Yathartha22

Yathartha22 Jul 4, 2018

Contributor

Is there any transcendental equation (other than something like cos(x) = 0) that can be solved/handled by transolve that is not an Add?

I can't think of any other than Add as of now if there are we can create another internal function for it.

It would be better to give a code snippet like

I will add the code snippet for better understanding.

Are different names used for the case when a solution is returned vs the case where a rewritten form is returned?

Yes, when the equation is recasted then new_eq will be updated and then solveset will be called to update the result, whereas when solution is returned the result is directly updated

invocation condition is placed inside `Add` case.
Once done, add a call to its corresponding *solving* *helper*.
For the class of equation define *identification* and

This comment has been minimized.

@smichr

smichr Jul 4, 2018

Member

The previous paragraph and this one need to be edited. The last one talks about th esolving helper but this paragraph sound like it is trying to introduce the helpers. I think there are three steps (if there are more than Add-type transcendentals):

  1. identify class
  2. write identifier
  3. write solvers

This comment has been minimized.

@Yathartha22

Yathartha22 Jul 4, 2018

Contributor

Okay explaining using three steps would be better:

Adding a new class of equation solver is a three-step procedure:

- Identify the type of the equations

Determine the type of the class of equations to which they belong, it could be of `Add`, `Pow`, 

etc. types. Seperate internal functions are used for each type, if it is already present include

identification and solving helpers within the function otherwise add a new internal function and then 

include identification and solving helpers within it.

### code snippet can be added here

- Define the identifying helpers

- Define the solving helpers

Something like this is okay?

lhs, rhs = list(ordered(f.args))
log_type_equation = expand_log(log(lhs)) - expand_log(log(-rhs))
return log_type_equation

This comment has been minimized.

@smichr

smichr Jul 4, 2018

Member

What is going to happen when the user selected the complex domain for the solution in solveset? The symbol won't be real, will it? And then the equation won't simplify as desired. I think you will get an infinite loop. But if there is no solution with numbers having an imaginary part then the real solution should be returned...and somewhere that change to a real variable will have to take place. I suspect that when solving in the complex domain an attempt at the complex solution will have to be followed with an attempt at a real solution.

This comment has been minimized.

@Yathartha22

Yathartha22 Jul 4, 2018

Contributor

What is going to happen when the user selected the complex domain for the solution in solveset?

The equation won't get simplified (by expand_log) but still it will be transformed to log hence it won't fall into infinite loop.

As of now untill _transolve gets finalised, @aktech suggested to not focus on solving in complex domain. I guess once this PR gets into a proper shape, I can add solving in complex domain in a seperate PR. What do you say?

This comment has been minimized.

@aktech

aktech Jul 9, 2018

Member

Not solving in complex domain is fine, but it should not break with an infinite loop.

result_rational = _solve_as_rational(equation, symbol, domain)
if isinstance(result_rational, ConditionSet):
# may be a transcendental type equation
result += _transolve(equation, symbol, domain)

This comment has been minimized.

@smichr

smichr Jul 4, 2018

Member

If you aren't going to implement solving in the complex domain then you should add the test here, perhaps:

...
# may be a transcendental type equation
if not domain.is_subset(S.Reals):
    raise NotImplementedError('solving transcendent equation in the complex domain.')
result += _transolve(...

This comment has been minimized.

@Yathartha22

Yathartha22 Jul 4, 2018

Contributor

Is getting a ConditionSet not a good idea?
As of now for the complex domain, ConditionSet is returned.

This comment has been minimized.

@smichr

smichr Jul 7, 2018

Member

oh, sorry. ConditionSet is fine.

@@ -978,6 +985,330 @@ def _solveset(f, symbol, domain, _check=False):
return result
def make_expr_args(f):

This comment has been minimized.

@Yathartha22

Yathartha22 Jul 7, 2018

Contributor

As of now the result of the routine will be used by identifying helpers.

This comment has been minimized.

@smichr

smichr Jul 7, 2018

Member

I suggest naming this term_factors since it returns the factors of all terms (contrary to what the current docstring says).

This comment has been minimized.

@Yathartha22

Yathartha22 Jul 8, 2018

Contributor

I don't think it returns factors of the expression, it just gives the terms as it is present in the expression

>>> make_expr_args(x**2 - 1)
[-1, x**2]

This comment has been minimized.

@Yathartha22

Yathartha22 Jul 8, 2018

Contributor

Also will it be a good idea to have this method placed somewhere else, this could be used as a public function (I don't think sympy has anything of this sought).

This comment has been minimized.

@smichr

smichr Jul 8, 2018

Member

don't think it returns factors of the expression

it doesn't return factors of expression, it returns factors of the terms of an expression:

>>> make_expr_args(v**x*y+z*(x+1))
[y, v**x, z, x + 1]

This comment has been minimized.

@Yathartha22

Yathartha22 Jul 8, 2018

Contributor

Oh sorry :), I miss understood a bit. I will update the name as well as the docstring

@@ -985,9 +985,10 @@ def _solveset(f, symbol, domain, _check=False):
return result
def make_expr_args(f):
def term_factors(f):

This comment has been minimized.

@smichr

smichr Jul 9, 2018

Member

What do you think about making this an iterator? So instead of returning all, yield one at a time. That could be used like uniq(term_factors(foo)) to give the unique term factors of an expression in an efficient manner.

This comment has been minimized.

@Yathartha22

Yathartha22 Jul 9, 2018

Contributor

Sounds good @smichr. I will make changes.

>>> from sympy import symbols
>>> from sympy.solvers.solveset import term_factors
>>> x = symbols('x')
>>> term_factors_obj = term_factors(x**2 - 1)

This comment has been minimized.

@smichr

smichr Jul 9, 2018

Member
>>> list(term_factors(-2 - x**2 + x*(x + 1)))
[-2, -1, x**2, x, x + 1]
>>> for term in term_factors_obj:
... print(term)
-1
x**2
"""
args = []

This comment has been minimized.

@smichr

smichr Jul 9, 2018

Member

not needed anymore

@aktech

For all the function that, you have created irrespective of private or public, add the following type of params and return type (for consistency):

    Parameters
    ==========

    f : Expr or a relational.
        The target equation or inequality
    symbol : Symbol
        The variable for which the equation is solved
    domain : Set
        The domain over which the equation is solved

    Returns
    =======

    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.

    Raises
    ======

    NotImplementedError
        The algorithms to solve inequalities in complex domain  are
        not yet implemented.
    ValueError
        The input is not valid.
    RuntimeError
        It is a bug, please report to the github issue tracker.

can be transformed to
`log(a) + f(x)*log(b) - log(c) - g(x)*log(d) = 0`
(under certain assumptions) and this can be solved with `solveset`
if `f(x)` and `g(x)` are polynomial in form.

This comment has been minimized.

@aktech

aktech Jul 9, 2018

Member

*are in polynomial form.

@@ -978,6 +985,338 @@ def _solveset(f, symbol, domain, _check=False):
return result
def term_factors(f):

This comment has been minimized.

@aktech

aktech Jul 9, 2018

Member

please make it private

lhs, rhs_s = invert_complex(f, 0, symbol, domain)
if isinstance(rhs_s, FiniteSet):
rhs = rhs_s.args[0]

This comment has been minimized.

@aktech

aktech Jul 9, 2018

Member

That's fine, but add an assert here to make your intentions clear, so that if someone changes that in solveset tomorrow, transolve should tell the world that it broke.

@@ -10,14 +10,17 @@
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, Equality,
Add, Mul)

This comment has been minimized.

@aktech

aktech Jul 9, 2018

Member

Mul imported twice

from sympy.core.containers import Tuple
from sympy.core.facts import InconsistentAssumptions
from sympy.core.numbers import I, Number, Rational, oo
from sympy.core.function import (Lambda, expand_complex, AppliedUndef, Function)
from sympy.core.function import (Lambda, expand_complex, AppliedUndef, Function,

This comment has been minimized.

@aktech

aktech Jul 9, 2018

Member

Please remove unused import Function

@aktech

This comment has been minimized.

Member

aktech commented Jul 9, 2018

Update the docstring at the top of the file solveset.py to include the solving of transcendental equations.

@XFAIL
def test_uselogcombine_2():

This comment has been minimized.

@aktech

aktech Jul 9, 2018

Member

Instead of 1 or 2 as function name suffix, write the difference between them as function name.

This comment has been minimized.

@Yathartha22

Yathartha22 Jul 10, 2018

Contributor

I don't see a difference here, these tests were already present I just moved them here.
I would rather include all in one.

@aktech

This comment has been minimized.

Member

aktech commented Jul 16, 2018

The PR looks good to me now. @smichr what do you say?

I will merge this in next 48 hours if no objections raised.

@smichr

This comment has been minimized.

Member

smichr commented Jul 16, 2018

The only question I had about design was whether there is any transcendental equation that is not of Add type that this solver will be attempting to solve? (I though I asked this but don't see it in the comments that are still showing.)

@Yathartha22

This comment has been minimized.

Contributor

Yathartha22 commented Jul 17, 2018

The only question I had about design was whether there is any transcendental equation that is not of Add type that this solver will be attempting to solve?

_transolve will be seeing the lhs of the equation of being Add, Mul, or Pow type. Currently, I have included only the add type because only solving exponential equations is implemented (which is basically of the add type).
There are Pow type of lhs like x**x - 2 = 0 (which can be solved in lambert form). So we need to include Pow case while implementing the lambert solver.
Similarly for Mul we have x*exp(x) - 1 = 0.
the thing is for this particular PR there doesn't seems a case, but needs to be taken care while implementing logarithmic and lambert solvers

@smichr

This comment has been minimized.

Member

smichr commented Jul 17, 2018

seeing the lhs of the equation

thanks for recap. +1

@aktech

aktech approved these changes Jul 17, 2018

Looks good to me.

@aktech aktech merged commit 7fe2ed3 into sympy:master Jul 17, 2018

1 check passed

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

This comment has been minimized.

Member

aktech commented Jul 17, 2018

@Yathartha22 thanks for the work! @smichr thanks for reviewing! This is in.

@Yathartha22

This comment has been minimized.

Contributor

Yathartha22 commented Jul 18, 2018

@aktech @smichr thanks for helping me out with this PR :)

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