nonlinsolve_till_24june_2016.py

##############################################################################
# ------------------------------nonlinsolve ---------------------------------#
##############################################################################


def substitution(system, symbols=None, result=[{}], known_symbols=[],
                 exclude=[], all_symbols=None):
    r""" Solves the `system` using substitution method.
    A helper function for `nonlinsolve`. This will be called from
    `nonlinsolve` when any equation(s) is non polynomial equation.

    Parameters
    ==========

    system : list of equations
        The target system of equations
    symbols : list of symbols to be solved.
        The variable(s) for which the system is solved
    known_symbols : list of solved symbols
        Values are known for these variable(s)
    result : An empty list or list of dict
        If No symbol values is known then empty list otherwise
        symbol as keys and corresponding value in dict.
    exclude : Set of expression.
        Mostly denominator expression(s) of the equations of the system.
        Final solution should not satisfy these expressions.
    all_symbols : known_symbols + symbols(unsolved).

    Returns
    =======

    A FiniteSet of ordered tuple of values of `all_symbols` for which the
    `system` has solution. Order of value is same as symbols order in the
    `all_symbols` list.

    Please note that general FiniteSet is unordered, the solution returned
    here is not simply a FiniteSet of solutions, rather it is a FiniteSet of
    ordered tuple, i.e. the first & only argument to FiniteSet is a tuple of
    solutions, which is ordered, & hence the returned solution is ordered.

    Also note that solution could also have been returned as an ordered tuple,
    FiniteSet is just a wrapper `{}` around the tuple. It has no other
    significance except for the fact it is just used to maintain a consistent
    output format throughout the solveset.

    Examples
    ========

    >>> from sympy.core.symbol import symbols
    >>> x, y = symbols('x, y', real = True)
    >>> from sympy.solvers.solveset import substitution
    >>> substitution([x + y], [x], [{y: 1}], [y], set([]), [x, y])
    {(-1, 1)}

    * when you want soln should not satisfy eq `x + 1 = 0`

    >>> substitution([x + y], [x], [{y: 1}], [y], set([x + 1]), [y, x])
    EmptySet()
    >>> substitution([x + y], [x], [{y: 1}], [y], set([x - 1]), [y, x])
    {(1, -1)}
    >>> substitution([x + y - 1, y - x**2 + 5], [x, y])
    {(-3, 4), (2, -1)}

    * Returns both real and complex solution

    >>> x, y, z = symbols('x, y, z')
    >>> from sympy import exp, sin
    >>> substitution([exp(x) - sin(y), y**2 - 4], [x, y])
    {(log(sin(2)), 2), (ImageSet(Lambda(_n, I*(2*_n*pi + pi) +
        log(sin(2))), Integers()), -2), (ImageSet(Lambda(_n, 2*_n*I*pi +
        Mod(log(sin(2)), 2*I*pi)), Integers()), 2)}

    >>> eqs = [z**2 + exp(2*x) - sin(y), -3 + exp(-y)]
    >>> substitution(eqs, [y, z])
    {(-log(3), -sqrt(-exp(2*x) - sin(log(3)))),
    (-log(3), sqrt(-exp(2*x) - sin(log(3)))),
    (ImageSet(Lambda(_n, 2*_n*I*pi + Mod(-log(3), 2*I*pi)), Integers()),
    ImageSet(Lambda(_n, -sqrt(-exp(2*x) + sin(2*_n*I*pi +
    Mod(-log(3), 2*I*pi)))), Integers())),
    (ImageSet(Lambda(_n, 2*_n*I*pi + Mod(-log(3), 2*I*pi)), Integers()),
    ImageSet(Lambda(_n, sqrt(-exp(2*x) + sin(2*_n*I*pi +
        Mod(-log(3), 2*I*pi)))), Integers()))}

    """

    from sympy.core.compatibility import ordered, default_sort_key
    from sympy import Complement

    if not system:
        return S.EmptySet

    if not symbols:
        raise ValueError('Symbols must be given, for which solution of the '
                         'system is to be found.')

    if hasattr(symbols[0], '__iter__'):
        symbols = symbols[0]

    try:
        sym = symbols[0].is_Symbol
    except AttributeError:
        sym = False

    if not sym:
        raise ValueError('Symbols or iterable of symbols must be given as '
                         'second argument, not type \
            %s: %s' % (type(symbols[0]), symbols[0]))

    # By default `all_symbols` will same as `symbols`
    if all_symbols is None:
        all_symbols = symbols

    def _unsolved_syms(eq, sort=False):
        """Returns the unsolved symbol present
        in the equation `eq`.
        """
        free = eq.free_symbols
        unsolved = (free - set(known_symbols)) & set(all_symbols)
        if sort:
            unsolved = list(unsolved)
            unsolved.sort(key=default_sort_key)
        return unsolved

    old_result = result
    # storing complements and intersection for particular symbol
    complements = {}
    intersections = {}

    def _solve_using_known_values(result, solver):
        """Solves the system using already known solution
        (result contains the dict <symbol: value>).
        solver is `solveset_complex` or `solveset_real`.
        """
        # stores imageset <expr: imageset(Lambda(n, expr), base)>.
        soln_imageset = {}
        # sort such that equation with the fewest potential symbols is first.
        # means eq with less variable first
        for index, eq in enumerate(
                ordered(system, lambda _: len(_unsolved_syms(_)))):
            newresult = []
            got_symbol = set()  # symbols solved in one iteration
            hit = False
            original_imageset = {}
            # if imageset expr is used to solve other symbol
            imgset_yes = False
            for res in result:
                if soln_imageset:
                    # find the imageset and use it's expr.
                    for key_res, value_res in res.items():
                        if isinstance(value_res, ImageSet):
                            res[key_res] = value_res.lamda.expr
                            original_imageset[key_res] = value_res
                            dummy_n = value_res.lamda.expr.atoms(Dummy).pop()
                            base = value_res.base_set
                            imgset_yes = (dummy_n, base)
                soln_imageset = {}
                # update eq with everything that is known so far
                eq2 = eq.subs(res)
                unsolved_syms = _unsolved_syms(eq2, sort=True)
                if not unsolved_syms:
                    if res:
                        newresult.append(res)
                    break  # skip as it's independent of desired symbols
                for sym in unsolved_syms:
                    not_solvable = False
                    try:
                        soln = solver(eq2, sym)
                        soln_new = S.EmptySet
                        if isinstance(soln, Complement):
                            # extract solution and complement
                            complements[sym] = soln.args[1]
                            soln = soln.args[0]
                            # complement will be added at the end
                        if isinstance(soln, Intersection):
                            # sometimes solveset returns Intersection with
                            # S.Reals when #11174 is fixed,then
                            # this if block should be removed

                            # Interval will be at 0th index
                            if soln.args[0] != Interval(-oo, oo):
                                intersections[sym] = soln.args[0]
                            soln_new += soln.args[1]
                        soln = soln_new if soln_new else soln
                        if index > 0 and solver == solveset_real:
                            # 1 symbol's real soln , another symbol may have
                            # corresponding complex soln.

                            if not isinstance(soln, ImageSet):
                                soln += solveset_complex(eq2, sym)
                    except NotImplementedError:
                        continue
                    if isinstance(soln, ConditionSet):
                            soln = S.EmptySet
                            # dont do `continue` we may get soln
                            # in terms of other symbol(s)
                            not_solvable = True
                    if isinstance(soln, Complement):
                            # extract solution and complement
                            complements[sym] = soln.args[1]
                            soln = soln.args[0]
                            # complement will be added at the end
                    if isinstance(soln, ImageSet):
                        soln_imagest = soln
                        expr2 = soln.lamda.expr
                        soln = FiniteSet(expr2)
                        soln_imageset[expr2] = soln_imagest

                    # if there is union of Imageset in soln
                    # no testcase for this if block
                    if isinstance(soln, Union):
                        soln_args = soln.args
                        soln = []
                        # need to use list.
                        # (finiteset union imageset).args is
                        # (finiteset, Imageset). But we want only iteration
                        # so append finteset elements and imageset.
                        for soln_arg2 in soln_args:
                            if isinstance(soln_arg2, FiniteSet):
                                list_finiteset = list(
                                    [sol_1 for sol_1 in soln_arg2])
                                soln += list_finiteset
                            else:
                                # ImageSet, Intersection, complement
                                # append directly
                                soln.append(soln_arg2)

                    for sol in soln:
                        # if there is union, then need to check
                        # complement, intersection, Imageset
                        # order should not be changed.
                        if isinstance(soln, Complement):
                            # extract solution and complement
                            complements[sym] = soln.args[1]
                            soln = soln.args[0]
                            # complement will be added at the end
                        if isinstance(sol, Intersection):
                            # sometimes solveset returns Intersection with
                            # S.Reals when #11174 is fixed,then
                            # this if block should be removed

                            # Interval will be at 0th index
                            if sol.args[0] != Interval(-oo, oo):
                                intersections[sym] = sol.args[0]
                            sol = sol.args[1]
                        # after intersection and complement Imageset should
                        # be checked.
                        if isinstance(sol, ImageSet):
                            soln_imagest = sol
                            sol = sol.lamda.expr
                            soln_imageset[sol] = soln_imagest

                        free = sol.free_symbols
                        if got_symbol and any([
                            ss in free for ss in got_symbol
                        ]):
                            # sol depends on previously solved symbols
                            # then continue
                            continue
                        rnew = res.copy()
                        # put each solution in res and append the new  result
                        # in the new result list (solution for symbol `s`)
                        # along with old results.
                        for k, v in res.items():
                            if isinstance(v, Expr):
                                # if any unsolved symbol is present
                                # Then subs known value
                                rnew[k] = v.subs(sym, sol)
                        # and add this new solution
                        rnew[sym] = sol
                        # for check solve use imageset expr
                        if not any(checksol(d, rnew) for d in exclude):
                            # if sol was imageset then add imageset
                            local_n = None
                            if imgset_yes:
                                local_n = imgset_yes[0]
                                base = imgset_yes[1]
                                # use ImageSet, we have dummy in sol
                                lam = Lambda(local_n, sol)
                                rnew[sym] = ImageSet(lam, base)
                            if soln_imageset:
                                rnew[sym] = soln_imageset[sol]
                                img_expr = soln_imageset[sol].lamda.expr
                                local_n = img_expr.atoms(Dummy).pop()
                            # restore original imageset
                            restore_sym = set(rnew.keys()) & \
                                set(original_imageset.keys())
                            dummy_n = None
                            base_img = None
                            for key_sym in restore_sym:
                                img = original_imageset[key_sym]
                                rnew[key_sym] = img
                                expr = img.lamda.expr
                                dummy_n = expr.atoms(Dummy).pop()
                                base_img = img.base_set
                            if local_n and restore_sym:
                                # also it means soln for `sym` contains
                                # dummy `n`, put this expr in ImageSet
                                # use only one dummy
                                img_expr = rnew[sym].lamda.expr
                                lam = Lambda(
                                    dummy_n,
                                    img_expr.subs(local_n, dummy_n)
                                )
                                rnew[sym] = ImageSet(lam, base_img)
                            # now append this solution
                            newresult.append(rnew)
                    hit = True
                    # solution got for s
                    if not not_solvable:
                        got_symbol.add(sym)
                if not hit:
                    raise NotImplementedError('could not solve %s' % eq2)
                else:
                    result = newresult
        return result
    # end def _solve_using_know_values

    new_result_real = _solve_using_known_values(old_result, solveset_real)
    new_result_complex = _solve_using_known_values(
        old_result, solveset_complex)

    # overall result
    result = new_result_real + new_result_complex
    result_finiteset = []
    for res in result:
        if not res:
            # if {None : None} is present.
            # No solution then first will be {None : None}
            continue
        # If length < len(all_symbols) means infinite soln.
        # Some or all the soln is dependent on 1 symbol.
        # eg. {x: y+2} then final soln is {x: y+2, y: y}
        if len(res) < len(all_symbols):
            solved_symbols = res.keys()
            unsolved = list(filter(
                lambda x: x not in solved_symbols, all_symbols))
            rcopy = res.copy()
            for unsolved_sym in unsolved:
                rcopy[unsolved_sym] = unsolved_sym
            # if symbol is not in `v` (unsolved symbol),
            # means it can take any value.
            # eg. if k, v => y, exp(x) then above lines will add {x: x}
            # if we have another symbol `z` then add {z: z} using above line.
            res = rcopy
        result_finiteset.append(res)

    if intersections:
        # If solveset have returned some intersection for any symbol
        # mostly intersection is Interval or Set
        result = []
        for res in result_finiteset:
            res_copy = res
            for key_res, value_res in res.items():
                for key_intersec, value_intersec in intersections.items():
                    if key_intersec == key_res:
                        res_copy[key_res] = \
                            Intersection(FiniteSet(value_res), value_intersec)
            result.append(res_copy)
        result_finiteset = result

    if complements:
        # If solveset have returned some complements for any symbol
        # Mostly complements are Set.
        result = []
        for res in result_finiteset:
            res_copy = res
            for key_res, value_res in res.items():
                for key_complement, value_complement in complements.items():
                    if key_complement == key_res:
                        res_copy[key_res] = \
                            FiniteSet(value_res) - value_complement
            result.append(res_copy)
        result_finiteset = result

    # convert to ordered tuple
    result = S.EmptySet
    for r in result_finiteset:
        temp = [r[sym] for sym in all_symbols]
        result += FiniteSet(tuple(temp))
    return result


def nonlinsolve(system, symbols):
    r"""
    Solve system of N non linear equations with M variables, which means both
    under and overdetermined systems are supported. Positive dimensional
    system is also supported (A system with infinitely many solutions is said
    to be positive-dimensional). In Positive dimensional system solution will
    be dependent on at least one symbol. Returns both real solution
    and complex solution. The possible number of solutions is zero,
    one or infinite.

    Parameters
    ==========

    system : list of equations
        The target system of equations
    symbols : list of Symbols
        symbols should be given as a sequence eg. list

    Returns
    =======

    A FiniteSet of ordered tuple of values of `symbols` for which the `system`
    has solution. Order of value is same as symbols order in the
    `all_symbols`list.

    Please note that general FiniteSet is unordered, the solution returned
    here is not simply a FiniteSet of solutions, rather it is a FiniteSet of
    ordered tuple, i.e. the first & only argument to FiniteSet is a tuple of
    solutions, which is ordered, & hence the returned solution is ordered.

    Also note that solution could also have been returned as an ordered tuple,
    FiniteSet is just a wrapper `{}` around the tuple. It has no other
    significance except for the fact it is just used to maintain a consistent
    output format throughout the solveset.

    For the given set of Equations, the respective input types
    are given below:

    .. math:: x*y - 1 = 0
    .. math:: 4*x**2 + y**2 - 5 = 0

    `system  = [x*y - 1, 4*x**2 + y**2 - 5]`
    `symbols = [x, y]`

    Examples
    ========

    >>> from sympy.core.symbol import symbols
    >>> from sympy.solvers.solveset import nonlinsolve
    >>> x, y, z = symbols('x, y, z', real=True)
    >>> nonlinsolve([x*y - 1, 4*x**2 + y**2 - 5], [x, y])
    {(-1, -1), (-1/2, -2), (1/2, 2), (1, 1)}

    * Positive dimensional system and complements:

    >>> from sympy import pprint
    >>> from sympy.polys.polytools import is_zero_dimensional
    >>> a, b, c, d = symbols('a, b, c, d', real=True)
    >>> eq1 =  a + b + c + d
    >>> eq2 = a*b + b*c + c*d + d*a
    >>> eq3 = a*b*c + b*c*d + c*d*a + d*a*b
    >>> eq4 = a*b*c*d -1
    >>> system = [eq1, eq2, eq3, eq4]
    >>> is_zero_dimensional(system)
    False
    >>> pprint(nonlinsolve(system, [a, b, c, d]), use_unicode=False)
      -1       1
    {(---, -d, -, {d} \ {0})}
       d       d
    >>> nonlinsolve([(x+y)**2 - 4, x + y - 2], [x, y])
    {(-y + 2, y)}

    * If some of the equations are non polynomial eq then `nonlinsolve`
    will call `substitution` function and returns real and complex solutions .

    >>> from sympy import exp, sin
    >>> nonlinsolve([exp(x) - sin(y), y**2 - 4], [x, y])
    {(log(sin(2)), 2), (ImageSet(Lambda(_n, I*(2*_n*pi + pi) +
        log(sin(2))), Integers()), -2), (ImageSet(Lambda(_n, 2*_n*I*pi +
        Mod(log(sin(2)), 2*I*pi)), Integers()), 2)}


    * If Non linear polynomial and zero dimensional system, then it returns
    both real and complex solutions using `solve_poly_system`.

    >>> from sympy import sqrt
    >>> nonlinsolve([x**2 - 2*y**2 -2, x*y - 2], [x, y])
    {(-2, -1), (2, 1), (-sqrt(2)*I, sqrt(2)*I), (sqrt(2)*I, -sqrt(2)*I)}


    * `nonlinsolve` can solve some linear( zero or positive dimensional)system
    solve this system (because it is using `groebner` function to get the
    groebner basis and then `substitution` function basis as the system).
    But it is not recommended to solve linear system using `nonlinsolve`,
    because `linsolve` is better for all kind of linear system.

    >>> nonlinsolve([x + 2*y -z - 3, x - y - 4*z + 9 , y + z - 4], [x, y, z])
    {(3*z - 5, -z + 4, z)}

    * system having polynomial equations:

    >>> e1 = sqrt(x**2 + y**2) - 10
    >>> e2 = sqrt(y**2 + (-x + 10)**2) - 3
    >>> nonlinsolve((e1, e2), (x, y))
    {(191/20, -3*sqrt(391)/20), (191/20, 3*sqrt(391)/20)}
    >>> nonlinsolve([x**2 + 2/y - 2, x + y - 3], [x, y])
    {(1, 2), (1 + sqrt(5), -sqrt(5) + 2), (-sqrt(5) + 1, 2 + sqrt(5))}
    >>> nonlinsolve([x**2 + 2/y - 2, x + y - 3], [y, x])
    {(2, 1), (2 + sqrt(5), -sqrt(5) + 1), (-sqrt(5) + 2, 1 + sqrt(5))}

    How nonlinsolve is better than old solver `_solve_system` :
    ===========================================================

    * A positive dimensional system solver : nonlinsolve can return
    solution for positive dimensional system. It finds the
    Groebner Basis of the positive dimensional system(calling it as
    basis) then we can start solving equation(having least number of
    variable first in the basis) using solveset and substituting that
    solved solutions into other equation(of basis) to get solution in
    terms of minimum variables. Here the important thing is how we
    are substituting the known values and in which equations.

    * Real and Complex both solutions : nonlinsolve returns solution both real
    and complex solution. If all the equations in the system are polynomial
    then using `solve_poly_system` both real and complex solution is returned.
    If all the equations in the system are not polynomial equation then solves
    Polynomial equations(If any) using `solve_poly_system` and goes to
    `substitution` method with this solution and non polynomial equation(s),
    to solve for unsolved variables. Here to solve for particular variable
    solveset_real and solveset_complex is used. For both real and complex
    solution function `_solve_using_know_values` is used inside `substitution`
    function.(`substitution` function will be called when there is non
    polynomial equation(s) is present). When solution is valid then add its
    general solution in the final result. There may the case when
    `solve_poly_system` will not give complex solution (only real solution)
    then that soln will be missing.If we want to get complex solution in
    general form then use `substitution` function.(it always returns all the
    solutions).


    * Complements and Intersection will be added if any : nonlinsolve maintains
    dict for complements and Intersections. If solveset find complements or
    Intersection with any Interval or set in during the execution of
    `substitution` function ,then complement and Intersection for that
    variable is added before returning final solution.

    """
    from sympy.solvers.solvers import _invert as _invert_solver
    from sympy.polys.polytools import is_zero_dimensional, groebner

    if not system:
        return S.EmptySet

    if not symbols:
        raise ValueError('Symbols must be given, for which solution of the '
                         'system is to be found.')

    if hasattr(symbols[0], '__iter__'):
        symbols = symbols[0]

    try:
        sym = symbols[0].is_Symbol
    except AttributeError:
        sym = False

    if not sym:
        raise ValueError('Symbols or iterable of symbols must be given as '
                         'second argument, not type \
            %s: %s' % (type(symbols[0]), symbols[0]))
    polys = []
    polys_expr = []
    nonpolys = []
    denominators = set()
    for eq in system:
        # Store denom expression if it contains symbol
        denominators.update(_simple_dens(eq, symbols))
        # TODO : solveset `_invert is not for multiple symbol.
        # using old solver's _invert for now
        independent, dependent = _invert_solver(eq, *symbols)
        eq = dependent - independent
        eq = eq.as_numer_denom()[0]
        # this will remove sqrt if symbols are under it
        poly = eq.as_poly(*symbols, extension=True)
        if poly is not None:
            polys.append(poly)
            polys_expr.append(poly.as_expr())
        else:
            nonpolys.append(eq)

    if len(symbols) == len(polys):
        # If all the equations in the system is poly
        if is_zero_dimensional(system):
            try:
                # try to solve 0 dimensional system using `solve_poly_system`
                result = solve_poly_system(polys, *symbols)
                # May be some extra soln is added because
                # we did `poly = g.as_poly(*symbols, extension=True)`
                # need to check that and remove if not a soln
                result_update = S.EmptySet
                for res in result:
                    dict_sym_value = dict(list(zip(symbols, res)))
                    if all(checksol(eq, dict_sym_value) for eq in system):
                        result_update += FiniteSet(res)
                return result_update
            except:
                # Right now it doesn't fail for any polynomial system of
                # equation if `solve_poly_system` fails then substitution
                # method will handle it. pass it
                pass

        # positive dimensional system
        # substitution method with groebner basis of the system
        basis = groebner(polys, symbols, polys=True)
        new_system = []
        for poly_eq in basis:
            new_system.append(poly_eq.as_expr())
        # solved_symbols = []
        result = [{}]
        result = substitution(
            new_system, symbols, result, [],
            denominators, symbols)

    else:
        # If alll the equations are not polynomial.
        # Use substitution method for nonpoly + polys eq
        result = substitution(
            polys_expr + nonpolys, symbols, exclude=denominators)

    if result:
        return result
    else:
        return S.EmptySet