GSOC 2014: Solvers

[Harsh Gupta]

Hi, I'm Harsh Gupta I will be GSOC applicant this year. I want to discuss the solvers idea given on the Idea's page.https://github.com/sympy/sympy/wiki/GSoC-2014-Ideas#solvers

Some of the aspect of the idea are discussed at https://groups.google.com/forum/?fromgroups=#!starred/sympy/8TM8cnuzkG8.


Aaron Meurer said:

I think with TransformationSet we can do quite a bit. That handles
sets like {f(x) | x in A}. I think what is missing is the basic set
builder {x | P(x)}, where P(x) is a boolean predicate.

Matthew Rocklin said:

> Real issue here - how to represent some solutions (e.g. sin(x)==0).

In sets we would represent this with
In [1]: imageset(k, pi*k, S.Integers)
Out[1]: {π⋅k | k ∊ ℤ}

Implementing a general set builder will be very useful, we will need to define basic set operations like union

and intersection on them. We might use a syntax like `BuildSet(expr, sym_in_inputset, cond)` 

In [1]: BuildSet(pi*k, (k, S.Integers), True).intersect(0, 10)

Out[1]: {pi, 2*pi, 3*pi}

Matthew Rocklin said:

It sounds like maybe solve should return a Set.

I think it will be necessary to return set if we want to represent the solution of expressions like `floor(x) - 5`. 

See https://code.google.com/p/sympy/issues/detail?id=3975.

One problem I see with returning sets is that it can break a lot of existing code.

As mentioned in the idea page we need to have a method to check if solve returns all the solutions. For polynomials or expressions which are 

solved by converting to polynomials we can compare the degree of the polynomial to the number of solutions returned.

Other method can be verifying the number of the solutions by using the derivative of the function. Say we are given a *continuous* and *differentiable*  function f(x) and it's derivative w.r.t x is g(x),

then if g(x) has n solutions then f(x) cannot have more than n + 1 solutions. So, if the solve returns n + 1 solutions for f(x) then we are guaranteed that we have found all the solutions.

For this we will need a discontinuity finder and we will also have to make sure that we have found all the solutions for g(x), i.e., the derivative of f(x), which can be done recursively or by using other methods.

A third way can be using numerical solver. Say we use solve to find the solution of f(x) for some interval [a, b] then we can also run the numerical solver for f(x) on [a, b] 

if the number of solutions by numerical solver and symbolic solve matches then we can be pretty sure that we have found out all the solutions of f(x) in [a, b]. 

--

Harsh

[Aaron Meurer]

On Tue, Jan 21, 2014 at 4:55 AM, Harsh Gupta <gupta....@gmail.com> wrote:
> Hi, I'm Harsh Gupta I will be GSOC applicant this year. I want to discuss
> the solvers idea given on the Idea's
> page.https://github.com/sympy/sympy/wiki/GSoC-2014-Ideas#solvers

Great to hear it. As noted on the ideas page, this one will require a
good deal of thought to be done in the application, so let's start
discussing.

>
> Some of the aspect of the idea are discussed at
> https://groups.google.com/forum/?fromgroups=#!starred/sympy/8TM8cnuzkG8.
>
>
> Aaron Meurer said:
>
>> I think with TransformationSet we can do quite a bit. That handles
>> sets like {f(x) | x in A}. I think what is missing is the basic set
>> builder {x | P(x)}, where P(x) is a boolean predicate.
>
>
>
> Matthew Rocklin said:
>
>>> > Real issue here - how to represent some solutions (e.g. sin(x)==0).
>>
>> In sets we would represent this with
>> In [1]: imageset(k, pi*k, S.Integers)
>> Out[1]: {π⋅k | k ∊ ℤ}
>
>
> Implementing a general set builder will be very useful, we will need to
> define basic set operations like union
> and intersection on them. We might use a syntax like `BuildSet(expr,
> sym_in_inputset, cond)`
>
> In [1]: BuildSet(pi*k, (k, S.Integers), True).intersect(0, 10)
> Out[1]: {pi, 2*pi, 3*pi}

So probably a good idea would be to look at the kinds of things that
solve might potentially want to return, and see where the deficiencies
are in the sets module. Things can get more complicated if you
consider higher dimensional sets.

>
>
> Matthew Rocklin said:
>
>> It sounds like maybe solve should return a Set.
>
>
> I think it will be necessary to return set if we want to represent the
> solution of expressions like `floor(x) - 5`.
> See https://code.google.com/p/sympy/issues/detail?id=3975.
> One problem I see with returning sets is that it can break a lot of existing
> code.

Yes, this is a problem. I think we should create custom objects that
extend sets that allow to keep the current interface, like sols[0] or
sols[0][x]. The current interface is of course quite a few interfaces,
which is part of what makes it such a mess, but I think the dict=True
way is the cleanest.

>
>
> As mentioned in the idea page we need to have a method to check if solve
> returns all the solutions. For polynomials or expressions which are
> solved by converting to polynomials we can compare the degree of the
> polynomial to the number of solutions returned.
>
> Other method can be verifying the number of the solutions by using the
> derivative of the function. Say we are given a *continuous* and
> *differentiable* function f(x) and it's derivative w.r.t x is g(x),
> then if g(x) has n solutions then f(x) cannot have more than n + 1
> solutions. So, if the solve returns n + 1 solutions for f(x) then we are
> guaranteed that we have found all the solutions.
> For this we will need a discontinuity finder and we will also have to make
> sure that we have found all the solutions for g(x), i.e., the derivative of
> f(x), which can be done recursively or by using other methods.

How do you determine the points of discontinuity of a function? In
general, you need to use solve (e.g., to find when you are dividing by
0). Hopefully you wouldn't try to recursively solve the same
expression, or you would get nowhere, except infinite recursion. But
it does need the same functionality, because if you have expr = p/q,
you need to be sure you know *all* the zeros of q to be sure where
expr is discontinuous.

>
> A third way can be using numerical solver. Say we use solve to find the
> solution of f(x) for some interval [a, b] then we can also run the numerical
> solver for f(x) on [a, b]
> if the number of solutions by numerical solver and symbolic solve matches
> then we can be pretty sure that we have found out all the solutions of f(x)
> in [a, b].

This has limitations too. Among them are:

- Numerical solvers make no such guarantee either. They just try to
find solutions using whatever method and return them when they do.

- Numerical solvers can be unstable for certain kinds of expressions.

- If you have symbolic parameters, this doesn't work so well. At best,
you can try a multidimensional numerical solver, but your chances of
finding all solutions in a higher dimensional space are even worse.

But this can definitely be used to weed out cases where we don't find
all the solutions. I was personally envisioning a more cautious
algorithm that only returns true when it is absolutely sure at each
stage, but I'm unclear just how far we can get with such a thing using
our current heuristics.

An audit of the current solve code might be in order. In particular,
I'd like to know:

1. what are the different "solvers"? (if we split solve into "hints"
like with dsolve, these would be the different hints), and
2. which are algorithmically complete (i.e., we know they will give
all solutions, or they can detect somehow if they may have missed
one)?

And this may raise auxiliary questions, like:

- to what degree can the different solvers be separated? For instance,
one solver (I'm not sure if it's actually implemented) would use
decompose() to solve recursively. How would such "recursive solvers"
look in a hints system?

- of those that are heuristic (not algorithmically complete), can they
be improved?

Another thing I'd like to know is if there's literature on solving
algorithms, particularly solving transcendental equations, and very
particularly on if there are any complete algorithms out there for
some class of equations.

By the way, one algorithm, relating to solving, which we don't have
implemented yet is CAS (cylindrical algebraic decomposition). This
would be a GSoC project in itself to implement, though.

Aaron Meurer

>
> --
> Harsh
>
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> For more options, visit https://groups.google.com/groups/opt_out.

[Matthew Rocklin]

Do we know how other computer algebra systems solve this problem?  How robust are the algorithms behind wolframalpha.com ?

[Harsh Gupta]

>> Great to hear it. As noted on the ideas page, this one will require a
>> good deal of thought to be done in the application, so let's start
>> discussing.

Thanks a lot, and sorry for the late reply

>> Another thing I'd like to know is if there's literature on solving
>> algorithms, particularly solving transcendental equations, and very
>> particularly on if there are any complete algorithms out there for
>> some class of equations.

I found a old paper called "SOLVING SYMBOLIC EQUATIONS WITH PRESS"
http://www.research.ed.ac.uk/portal/files/413486/Solving_Symbolic_Equations_%20with_PRESS.pdf

>> Do we know how other computer algebra systems solve this problem? How robust are the algorithms behind wolframalpha.com ?

I have found another paper "A Review of Symbolic Solvers"
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.44.9444&rep=rep1&type=pdf
and according to it Mathematica performs performs pretty bad.

>> An audit of the current solve code might be in order. In particular,
>> I'd like to know:
>>
>> 1. what are the different "solvers"? (if we split solve into "hints"
>> like with dsolve, these would be the different hints), and
>> 2. which are algorithmically complete (i.e., we know they will give
>> all solutions, or they can detect somehow if they may have missed
>> one)?
>>
>> And this may raise auxiliary questions, like:
>>
>> - to what degree can the different solvers be separated? For instance,
>> one solver (I'm not sure if it's actually implemented) would use
>> decompose() to solve recursively. How would such "recursive solvers"
>> look in a hints system?
>>
>> - of those that are heuristic (not algorithmically complete), can they
>> be improved?

I'm going through the solvers code and will answer these questions soon.

[Aaron Meurer]

On Fri, Jan 24, 2014 at 2:02 PM, Harsh Gupta <gupta....@gmail.com> wrote:
>>> Great to hear it. As noted on the ideas page, this one will require a
>>> good deal of thought to be done in the application, so let's start
>>> discussing.
>
> Thanks a lot, and sorry for the late reply
>
>>> Another thing I'd like to know is if there's literature on solving
>>> algorithms, particularly solving transcendental equations, and very
>>> particularly on if there are any complete algorithms out there for
>>> some class of equations.
>
> I found a old paper called "SOLVING SYMBOLIC EQUATIONS WITH PRESS"
> http://www.research.ed.ac.uk/portal/files/413486/Solving_Symbolic_Equations_%20with_PRESS.pdf
>
>>> Do we know how other computer algebra systems solve this problem? How robust are the algorithms behind wolframalpha.com ?
>
> I have found another paper "A Review of Symbolic Solvers"
> http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.44.9444&rep=rep1&type=pdf
> and according to it Mathematica performs performs pretty bad.

That was in 1996.

Nonetheless this, along with the Wester paper, should provide some
good test cases so we can see what can be done that we can't do.

Aaron Meurer

>
>>> An audit of the current solve code might be in order. In particular,
>>> I'd like to know:
>>>
>>> 1. what are the different "solvers"? (if we split solve into "hints"
>>> like with dsolve, these would be the different hints), and
>>> 2. which are algorithmically complete (i.e., we know they will give
>>> all solutions, or they can detect somehow if they may have missed
>>> one)?
>>>
>>> And this may raise auxiliary questions, like:
>>>
>>> - to what degree can the different solvers be separated? For instance,
>>> one solver (I'm not sure if it's actually implemented) would use
>>> decompose() to solve recursively. How would such "recursive solvers"
>>> look in a hints system?
>>>
>>> - of those that are heuristic (not algorithmically complete), can they
>>> be improved?
>
> I'm going through the solvers code and will answer these questions soon.
>

[Harsh Gupta]

I'm reading and understanding the solvers code. I have started
documenting it here https://github.com/sympy/sympy/wiki/solvers.

@Matthew
For implementing and dealing with infinite sets I've found a draft by
Richard Fateman
http://www.cs.berkeley.edu/~fateman/papers/sets.pdf

I have skimmed through it and it appears all of the techniques
described there are implementable in sympy.

--
Harsh

[Aaron Meurer]

On Sun, Jan 26, 2014 at 8:11 PM, Harsh Gupta <gupta....@gmail.com> wrote:
> I'm reading and understanding the solvers code. I have started
> documenting it here https://github.com/sympy/sympy/wiki/solvers.

That's great. Once you are finished we should clean it up and convert
it to a page in our docs.

Aaron Meurer

[mario]

You wrote "Methods to solve solvable quintics are implemented in sympy."

Apparently only for some of them: it does not solve
``x**5 - 5*x**4 + 30*x**3 - 50*x**2 + 55*x - 21 = 0``

taken from http://en.wikipedia.org/wiki/Quintic_function

[Harsh Gupta]

> Apparently only for some of them: it does not solve
> ``x**5 - 5*x**4 + 30*x**3 - 50*x**2 + 55*x - 21 = 0``

Thanks. Yes, not all of them, Only equations of form x**5 + p*x**3 +
q*x**2 + r*x + s, no fourth order terms are solvable.
The implementation was added in
https://github.com/sympy/sympy/pull/1746. So, there is scope of
improvement. I wonder
how many of other methods of solving solvable quintics can be
implemented without a knowledge of abstract algebra.
Aaron Meurer can you guide me on this?

[mario]

One can always eliminate the fourth order by a translation, in this case by one; one obtains then
``x**5 + 20*x**3 + 20*x**2 + 30*x + 10``, which is solved by SymPy; it would be useful to automatize this step.

However the real solution is not simplified, it has ``count_ops = 267``, vs ``count_ops=11`` for the simplified
real solution ``S('2^(1/5) - 4^(1/5) + 8^(1/5) - 16^(1/5)')``

[Aaron Meurer]

Ah, apparently the paper
(http://www.ams.org/journals/mcom/1991-57-195/S0025-5718-1991-1079014-X/S0025-5718-1991-1079014-X.pdf
("Solving solvable quintics", D. S. Dummit)), glosses over this,
because of this fact. What is the algorithm to eliminate the fourth
order term?

How can we simplify the roots. Will the more advanced sqrtdenest
algorithms help (https://code.google.com/p/sympy/issues/detail?id=93)?

Aaron Meurer

[Aaron Meurer]

Oh I guess it's easy:

In [32]: print(collect(a.subs(x, x - y).expand(), x))
x**5 + x**4*(-5*y - 5) + x**3*(10*y**2 + 20*y + 30) + x**2*(-10*y**3 -
30*y**2 - 90*y - 50) + x*(5*y**4 + 20*y**3 + 90*y**2 + 100*y + 55) -
y**5 - 5*y**4 - 30*y**3 - 50*y**2 - 55*y - 21

So in this case, we just need to shift by y=-1.

Aaron Meurer

[mario]

I do not know how to simplify these roots.

[someone]

Hi,

> > How can we simplify the roots. Will the more advanced sqrtdenest
> > algorithms help
> > (https://code.google.com/p/sympy/issues/detail?id=93)?

I'm not sure if this is of any help here, but rechecking and
implementing the complete root denesting algorithm should be
a good thing in any case. See the work by R. Zippel and S. Landau.

I don't the state of the code in sympy.

[Aaron Meurer]

If the algorithm is actually complete, and works for higher order
roots (not just square roots), then I suppose it should work here.

Aaron Meurer

[Harsh Gupta]

Sorry for being passive for some days. I was busy with assignments and
some other stuff.

> See the work by R. Zippel and S. Landau.

Are you talking these articles:

http://www.sciencedirect.com/science/article/pii/074771719290004N ( A
note on “Zippel Denesting” by Susan Landau )
http://www.computer.org/csdl/proceedings/focs/1989/1982/00/063496.pdf
(Simplification of Nested Radicals by Susan Landau)
http://www.sciencedirect.com/science/article/pii/S0747717185800146 (
Simplification of expressions involving radicals by Richard Zippel)

I also found this article:
http://pdf.aminer.org/001/059/825/simplification_of_nested_radicals.pdf
(How to denest Ramanujan's Nested Radicals by Johannes Blomer)

I think having general denesting techniques will be very useful when
while working with higher order polynomials, we should try to
implement them in sympy. I read the first article (A note on Zipple
Denesting) and skimmed over the rest of them. It appears that all of
them need at least some knowledge of group theory in general. I have
no previous experience in working with group theory. So I want to know
what is the status of group theory and field theory in Sympy? Also I
wanted to know if it would be realistic if I read Group theory before
the summer and add it to the GSOC idea.

[Aaron Meurer]

You'll need an understanding of basic field theory to understand those
algorithms (maybe group theory too if they are using Galois theory).

Aaron Meurer

[Harsh Gupta]

As part of handling infinite solutions for something like (sin(x) ==
0) I have started some work in
https://github.com/sympy/sympy/pull/2904.
After this we will be able to have intersection of many infinite
imagesets with finite sets.
I think we should list out other set theoretic capabilities we will be
needing to handle infinite solutions.

BTW Matthew will you be a mentor for this project? The set theoretic
capabilities of sympy is mainly due to you and we will needing more of
them for this project. Also your Set Expression module will be useful
if want return a Set as a solution.

--
Harsh

[Matthew Rocklin]

Neither qualified nor particularly well qualified to mentor this project.  For me the sets module was a good example of how to cleanly represent and simplify expressions; I wasn't interested in the mathematical details of how such simplifications are accomplished.  I'm more interested in infrastructure than in actual mathematics.

That being said, I depend on SymPy's solvers quite a bit, and so strongly encourage projects that improve them.

[Aaron Meurer]

(This applies to all potential GSoC applicants)

Don't worry if there is someone to mentor your GSoC project. Just
apply, and if we like the application, we will find someone to mentor
it. We try very hard to avoid not accepting a good student just
because there are no mentors for the project.

Aaron Meurer

[someone]

Hi,

> >> >> Are you talking these articles:
> >> >>
> >> >> http://www.sciencedirect.com/science/article/pii/074771719290004N
> >> >> ( A note on "Zippel Denesting" by Susan Landau )
> >> >> http://www.computer.org/csdl/proceedings/focs/1989/1982/00/063496.pdf
> >> >> (Simplification of Nested Radicals by Susan Landau)
> >> >> http://www.sciencedirect.com/science/article/pii/S0747717185800146
> >> >> ( Simplification of expressions involving radicals by Richard
> >> >> Zippel)
> >> >>
> >> >> I also found this article:
> >> >> http://pdf.aminer.org/001/059/825/simplification_of_nested_radicals.pdf
> >> >> (How to denest Ramanujan's Nested Radicals by Johannes Blomer)

Yes. One should implement that. Maybe looking at the algorithm in the last
paper first because it seems to be the simplest one.

[Jigar Mistry]

i am also interested in solvers module. how can i contribute to simpy using this solvers module.Plz help me....

[Harsh Gupta]

I realize that there is less than a month left for the application
deadline and it has already been a month since I posted this thread.
So, I need to come up with design details fast. I have written some
details about using sets as output for solvers in this PR
https://github.com/sympy/sympy/pull/2948/files. PR because I feel it
is better platform to facilitate this discussion than the mailing
list. Please have a look.

[Harsh Gupta]

Hi, I have put my the first draft of the proposal on the sympy wiki.

https://github.com/sympy/sympy/wiki/GSoC-2014-Application-Harsh-Gupta:-Solvers

Please have a look, your suggestions will be valuable.

--
Harsh

[Aaron Meurer]

I think it looks good so far. I like how you cleanly organized your
pull requests. That helps a lot.

My comments:

- I'm putting this up front because I think it's the most important
point. I feel that you are only addressing half of the issue at
https://github.com/sympy/sympy/issues/6659, namely the output format
of solve. What about the input API? Here is the parameter list for
solve() (I had to construct this manually from the docstring since the
code uses **kwargs):

solve(f, *symbols, dict=False, set=False, exclude=(), check=True,
numerical=True, minimal=False, warning=False, simplify=True,
force=False, rational=True, manual=False, implicit=False,
minimal=False, quick=False)

That's assuming there aren't any undocumented parameters.

And don't get me started on the mess of linear equation solving
functions. We probably have a dozen functions that solve linear
equations or systems of linear equations (see some of the discussion
at https://github.com/sympy/sympy/pull/2580).

I like that you're adding new stuff and cleaning up the output API,
but you should think about how to clean up the input API as well.

- Your example

In[]: soln = solve(sin(x)*sin(y), (x, y), input_set = Interval(0,
4)*Interval(0, 4))

is a bit confusing to me. The input_set argument gives a 2-dimensional
set, but how are you to know which axis is x and which is y?

- You talk a lot about using sets, which I think is a good idea. But
you should think about how you can also use the assumptions. Maybe
there is a clean way that we can go back and forth between assumptions
and sets that requires minimal code duplication, and also allows each
to take advantage of the algorithms implemented in the other (by the
way, when I say assumptions, you should probably only worry about the
new assumptions, i.e., the stuff in the Q object).

- Your section on singularities doesn't really make it clear why we
need the solvers to be improved to do this, namely, that we need to
know exactly when we have all the solutions to not get wrong results.

- How will we handle that situation (finding all solutions)? What if
we can't say anything? Can we still represent objects in such a way
that it is not wrong (basically by somehow saying, "here are 'some' of
the solutions, but maybe not all of them", and ditto for anything that
uses solve, like singularities)? Maybe Piecewise is sufficient
somehow?

- I'm glad to see you are working on improving the simplicity of the
results of solve, such as radical denesting. A potential
simplification routine uses minpoly() on an algebraic expression and
then uses solve to find a simpler solution, but this relies on the
results of solve being simple. Do the radical denesting algorithms
work with symbolic entries as well?

- Did you plan to add any new solvers? I think there are still quite a
few cases that we can't solve. Some higher degree irreducible
polynomials for instance (not all higher degree polynomials are
solvable by radicals but some are). There will also be a lot to
implement once we are able to even represent the solutions to
sin(x)=0.

Aaron Meurer

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CADN8iuohUoie3-eQHG5C_%2BEtKsfs%2BJ%2Bk2P-_DdLMPSBbLtXXDA%40mail.gmail.com.
> For more options, visit https://groups.google.com/d/optout.

[Matthew Rocklin]

I think that it is far more important to build a sane and extensible framework around solve than to focus on building new solvers.  My guess is that this also requires a different skill set.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6KJVcvAhtX0iXPN%3DRy9DVrt642KmXLkAxLSxMDpNy5jiA%40mail.gmail.com.

[Aaron Meurer]

I agree, although it is usually a good idea to have a lot of
motivating examples to write a good framework. But in this case, the
existing solvers should be enough (with perhaps the exception of the
solvers that cannot exist currently because the current framework
doesn't allow them, like parameterized solutions).

Aaron Meurer

> https://groups.google.com/d/msgid/sympy/CAJ8oX-F3Q6psXPJLk3-ZnWJtq9F0z42fo5Dp4epBoX%2BUg9Zt7Q%40mail.gmail.com.

[Harsh Gupta]

> - Your example

> 

> In[]: soln = solve(sin(x)*sin(y), (x, y), input_set = Interval(0,

> 4)*Interval(0, 4))

> 

> is a bit confusing to me. The input_set argument gives a 2-dimensional

> set, but how are you to know which axis is x and which is y?

The axis is determined by the order of variables in which they appear in the

argument of solve. Defining the input set in this way will give us special

advantage of being able to take values from the sets which are traditionally

hard to define. For example if we want the variable to come from a circle we

can do it like.

In[]: solve(f(x, y), (x, y), input_set = imageset((x, y), x**2 + y**2 < 1, S.Reals*S.Reals))

For a L shape domain we can do

In[]: solve(f(x, y), (x, y), input_set = Intersection(Interval(0, 2)*Interval(0, 3), Interval(1, 2)*Interval(1, 3))

I cannot be sure that sets will be able to seamlessly handle such sets, but

I really think this API will scale.

> What about the input API?

>

> solve(f, *symbols, dict=False, set=False, exclude=(), check=True,

> numerical=True, minimal=False, warning=False, simplify=True,

> force=False, rational=True, manual=False, implicit=False,

> minimal=False, quick=False)

I think most of the flags are not needed. The flags like `dict` and `set`, won't be need as we

are unifying the output to set. Do we have any estimate of how many of these flags are actually used by users?

> - You talk a lot about using sets, which I think is a good idea. But

> you should think about how you can also use the assumptions. Maybe

> there is a clean way that we can go back and forth between assumptions

> and sets that requires minimal code duplication, and also allows each

> to take advantage of the algorithms implemented in the other (by the

> way, when I say assumptions, you should probably only worry about the

> new assumptions, i.e., the stuff in the Q object).

As mentioned in the comment by Matthew https://github.com/sympy/sympy/pull/2948#issuecomment-36592347

Assumptions can answer questions like if k is in N, is k*(k+1)/2 in N? This will

clearly help in resolving some of the set operations. btw I think assumptions

is wrong in this result.

In [20]: with assuming(Q.integer(k/2)):

    ...:     print(k in S.Integers)

    ...:     

False

 

> - How will we handle that situation (finding all solutions)? What if

> we can't say anything? Can we still represent objects in such a way

> that it is not wrong (basically by somehow saying, "here are 'some' of

> the solutions, but maybe not all of them", and ditto for anything that

> uses solve, like singularities)? Maybe Piecewise is sufficient

> somehow?

 

We will return a set as a solution if and only if we have found all the

solutions in the given domain. For every other case we will be using the

unevaluated Solve object. We will have a attribute in the unevaluated solve

named "know_solutions" to say "here are some solutions". Yes Picewise might be

helpful, but for now I can't think of a clear way to say "assuming this parameter 'a' is positive

the solutions are ..."

>  Do the radical denesting algorithms

> work with symbolic entries as well?

I don't think so.

> - Did you plan to add any new solvers? I think there are still quite a

> few cases that we can't solve. Some higher degree irreducible

> polynomials for instance (not all higher degree polynomials are

> solvable by radicals but some are). There will also be a lot to

> implement once we are able to even represent the solutions to

> sin(x)=0.

There was one algorithm

I discussed on the discussion https://github.com/sympy/sympy/pull/2948#issuecomment-36970134. 

By that algo we will be able to solve some special cases like

`sin(x) == x`. I will implement that and we might come up with new algorithm in the

process of rewriting current solvers. If time permits I'll surely try to

implement new solvers.

[Matthew Rocklin]

> input_set = imageset((x, y), x**2 + y**2 < 1, S.Reals*S.Reals))

Note that sets doesn't work this way, we don't take booleans.  I'm actually pretty sure you know this already, I just want to show off how sets works to others.

Instead you could do the following

In [4]: unit_disk = imageset(Lambda((r, theta), (r * cos(theta), r * sin(theta))), Interval(0, 1) * Interval(0, 2*pi))

In [5]: unit_disk

Out[5]: {(r⋅cos(θ), r⋅sin(θ)) | r, θ ∊ [0, 1] × [0, 2⋅π]}

In [6]: (0.5, 0.5) in unit_disk

Out[6]: True

In [7]: (0.5, 2) in unit_disk

Out[7]: False

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CADN8iuq9aVesJzFFE0wpG9_eJd6C7RW-7Lb%2BRZ3OHA4AXP1yXA%40mail.gmail.com.

[Aaron Meurer]

Quite a few people use dict=True, as that's the recommended format to
get a consistent result. As with any public API, we'd have to
deprecate it before removing it.

What about the other flags? And what about the dozen linear solvers? I
don't expect you to have the right answers now these things this (I'm
not even sure myself what should be done about them), but you should
be thinking about them.

>
>> - You talk a lot about using sets, which I think is a good idea. But
>> you should think about how you can also use the assumptions. Maybe
>> there is a clean way that we can go back and forth between assumptions
>> and sets that requires minimal code duplication, and also allows each
>> to take advantage of the algorithms implemented in the other (by the
>> way, when I say assumptions, you should probably only worry about the
>> new assumptions, i.e., the stuff in the Q object).
>
> As mentioned in the comment by Matthew
> https://github.com/sympy/sympy/pull/2948#issuecomment-36592347
> Assumptions can answer questions like if k is in N, is k*(k+1)/2 in N? This
> will
> clearly help in resolving some of the set operations. btw I think
> assumptions
> is wrong in this result.
>
> In [20]: with assuming(Q.integer(k/2)):
> ...: print(k in S.Integers)
> ...:
> False

The problem is that we always return an answer with __contains__. This
has nothing to do with the assumptions, but rather the sets module.
__contains__ should raise an exception if it can't determine (note
that Python forces "in" to always return a boolean, so we can't return
anything symbolic).

But beyond that, I don't think the sets use the new assumptions at all
(correct me if I am wrong).

I also noticed this:

In [32]: S.Reals.contains(x)
Out[32]: True

(x is just Symbol('x')).

>
>> - How will we handle that situation (finding all solutions)? What if
>> we can't say anything? Can we still represent objects in such a way
>> that it is not wrong (basically by somehow saying, "here are 'some' of
>> the solutions, but maybe not all of them", and ditto for anything that
>> uses solve, like singularities)? Maybe Piecewise is sufficient
>> somehow?
>
> We will return a set as a solution if and only if we have found all the
> solutions in the given domain. For every other case we will be using the
> unevaluated Solve object. We will have a attribute in the unevaluated solve
> named "know_solutions" to say "here are some solutions". Yes Picewise might
> be
> helpful, but for now I can't think of a clear way to say "assuming this
> parameter 'a' is positive
> the solutions are ..."

I guess one idea would be to union the set with some "additional
solutions" set, for which not much is known about. So the result, in
the case where we don't know if we have all the solutions, would be

SolutionSet U NotFoundSolutions

Then things that use the solutions like discontinuities would
propagate the union (discontinuities should return a Set object too
btw).

The unknown solutions set may have properties known, even if its exact
cardinality isn't. For instance, the zero set of a continuous function
is always closed.

>
>
>> Do the radical denesting algorithms
>> work with symbolic entries as well?
>
> I don't think so.

So in that case, one should try to improve the solvers themselves to
return simpler answers in the first place, if possible.

>
>
>> - Did you plan to add any new solvers? I think there are still quite a
>> few cases that we can't solve. Some higher degree irreducible
>> polynomials for instance (not all higher degree polynomials are
>> solvable by radicals but some are). There will also be a lot to
>> implement once we are able to even represent the solutions to
>> sin(x)=0.
>
> There was one algorithm
> I discussed on the discussion
> https://github.com/sympy/sympy/pull/2948#issuecomment-36970134.
> By that algo we will be able to solve some special cases like
> `sin(x) == x`. I will implement that and we might come up with new algorithm
> in the
> process of rewriting current solvers. If time permits I'll surely try to
> implement new solvers.

As I noted before, algorithms that just extend what is already there
should go at the end of the timeline. But algorithms that will provide
motivating examples for the solve API should be implemented sooner, at
least in part.

(p.s. don't forget to be updating your proposal with ideas from this discussion)

Aaron Meurer

>
>
> --
> You received this message because you are subscribed to the Google Groups
> "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to sympy+un...@googlegroups.com.
> To post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit

> https://groups.google.com/d/msgid/sympy/CADN8iuq9aVesJzFFE0wpG9_eJd6C7RW-7Lb%2BRZ3OHA4AXP1yXA%40mail.gmail.com.

[Matthew Rocklin]

S.Integers._contains(self, other) just calls ask(Q.integer(other))

S.Reals is actually just Interval(-oo, oo).  Interval._contains depends on the < and > operators.  It returns true because the following return true

In [9]: x < oo

Out[9]: True

In [10]: x > -oo

Out[10]: True

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6JtpThrtt7TVqJ162gNO0qepP_hf%3DKcARRaGdCEST%3Dkww%40mail.gmail.com.

[Aaron Meurer]

On Sun, Mar 16, 2014 at 2:27 PM, Matthew Rocklin <mroc...@gmail.com> wrote:
> S.Integers._contains(self, other) just calls ask(Q.integer(other))

It should raise an exception when ask() returns None.

>
> S.Reals is actually just Interval(-oo, oo). Interval._contains depends on
> the < and > operators. It returns true because the following return true
> In [9]: x < oo
> Out[9]: True
>
> In [10]: x > -oo
> Out[10]: True

It would be better to use ask and Q.real. The inequalities let you put
complex variables in them, because otherwise it would be annoying to
have to always use x = Symbol('x', real=True). But it leads to
incorrect results in this case.

By the way, I didn't know that the sets use exclusively the new
assumptions (assuming that is indeed what you are saying here). Given
what I've seen in the new assumptions, I'm not so sure this is such a
hot idea (namely, the handlers have a ton of mathematical
inaccuracies), but it's nice that they actually are used anywhere
nonetheless.

Aaron Meurer

> https://groups.google.com/d/msgid/sympy/CAJ8oX-E9B1uFoqw2_FTRxkJcm7jNk%3DM8NdKkco8bzsx3tYae5A%40mail.gmail.com.

[Harsh Gupta]

Here's what I think we should we with the parameters. The list is generated from the docstring.

'dict'=True (default is False)

        return list (perhaps empty) of solution mappings

'set'=True (default is False)

    return list of symbols and set of tuple(s) of solution(s)

Not needed, sets as output will take care of it.

'implicit=True (default is False)'

    allows solve to return a solution for a pattern in terms of

    other functions that contain that pattern; this is only

    needed if the pattern is inside of some invertible function

    like cos, exp, ....

We won't be able to support such output with our new output API

assert solve(x - exp(x), x, implicit=True) == [exp(x)]

'force=True (default is False)'

    make positive all symbols without assumptions regarding sign.

The user can define this in the input_set parameter.

'rational=True (default)'

    recast Floats as Rational; if this option is not used, the

    system containing floats may fail to solve because of issues

    with polys. If rational=None, Floats will be recast as

    rationals but the answer will be recast as Floats. If the

    flag is False then nothing will be done to the Floats.

I don't know if people use it. If such flag is not popular we will use

it. Anyway sympy solvers are mainly for "exact" symbolic coefficient so floating point

coefficients doesn't really fit in.

'exclude=[] (default)'

    don't try to solve for any of the free symbols in exclude;

    if expressions are given, the free symbols in them will

    be extracted automatically.

The new API makes is mandatory to enter the variable as an argument so all symbols

not given in the variables list is excluded.

'particular=True (default is False)'

    instructs solve to try to find a particular solution to a linear

    system with as many zeros as possible; this is very expensive

I'm not sure about this one.

'check=True (default)'

    If False, don't do any testing of solutions. This can be

    useful if one wants to include solutions that make any

    denominator zero.

'numerical=True (default)'

    do a fast numerical check if ``f`` has only one symbol.

'minimal=True (default is False)'

    a very fast, minimal testing.

'quick=True (default is False)'

    when using particular=True, use a fast heuristic instead to find a

    solution with many zeros (instead of using the very slow method

    guaranteed to find the largest number of zeros possible)

'warning=True (default is False)'

    show a warning if checksol() could not conclude.

'simplify=True (default)'

    simplify all but cubic and quartic solutions before

    returning them and (if check is not False) use the

    general simplify function on the solutions and the

    expression obtained when they are substituted into the

    function which should be zero

'manual=True (default is False)'

    do not use the polys/matrix method to solve a system of

    equations, solve them one at a time as you might "manually".

All these parameters trade the speed of the solvers with the accuracy. I know there have been issue regarding speed, see https://github.com/sympy/sympy/issues/6611. But I still think we should not be using these flags for the reason of maintaining simplicity. If a lot of people report that our solvers are too slow then we might reintroduce a few of them.

It turns out that we might go about not using any of the current parameters and have a clean simple input api as 

solve(<expression>, <variables>, input_set=<domain of variables>)

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6Kw671B%2B4XKcyV%2BeMLY3gffy%3DGqiAKrqVwjOEswXnRS6A%40mail.gmail.com.

For more options, visit https://groups.google.com/d/optout.

--
Harsh

[Aaron Meurer]

On Sun, Mar 16, 2014 at 3:20 PM, Harsh Gupta <gupta....@gmail.com> wrote:
> Here's what I think we should we with the parameters. The list is generated
> from the docstring.
>
> 'dict'=True (default is False)
> return list (perhaps empty) of solution mappings
>
> 'set'=True (default is False)
> return list of symbols and set of tuple(s) of solution(s)
>
> Not needed, sets as output will take care of it.
>
> 'implicit=True (default is False)'
> allows solve to return a solution for a pattern in terms of
> other functions that contain that pattern; this is only
> needed if the pattern is inside of some invertible function
> like cos, exp, ....
>
> We won't be able to support such output with our new output API
> assert solve(x - exp(x), x, implicit=True) == [exp(x)]

How is this implemented?

>
>
> 'force=True (default is False)'
> make positive all symbols without assumptions regarding sign.
>
> The user can define this in the input_set parameter.
>
> 'rational=True (default)'
> recast Floats as Rational; if this option is not used, the
> system containing floats may fail to solve because of issues
> with polys. If rational=None, Floats will be recast as
> rationals but the answer will be recast as Floats. If the
> flag is False then nothing will be done to the Floats.
>
> I don't know if people use it. If such flag is not popular we will use
> it. Anyway sympy solvers are mainly for "exact" symbolic coefficient so
> floating point
> coefficients doesn't really fit in.

The issue here is with the way the polys handles floating point
numbers. Basically, some things do not work if you use floating point
numbers but they do if you use rationals. Once the base issue is
fixed, the need for this flag will be obviated.

>
> 'exclude=[] (default)'
> don't try to solve for any of the free symbols in exclude;
> if expressions are given, the free symbols in them will
> be extracted automatically.
>
> The new API makes is mandatory to enter the variable as an argument so all
> symbols
> not given in the variables list is excluded.
>
> 'particular=True (default is False)'
> instructs solve to try to find a particular solution to a linear
> system with as many zeros as possible; this is very expensive
>
> I'm not sure about this one.

This is part of the messy linear system API that I alluded to. I'm
not clear what to do this this without seeing the full picture of all
the linear solvers.

>
> 'check=True (default)'
> If False, don't do any testing of solutions. This can be
> useful if one wants to include solutions that make any
> denominator zero.
>
> 'numerical=True (default)'
> do a fast numerical check if ``f`` has only one symbol.
>
> 'minimal=True (default is False)'
> a very fast, minimal testing.
>
> 'quick=True (default is False)'
> when using particular=True, use a fast heuristic instead to find a
> solution with many zeros (instead of using the very slow method
> guaranteed to find the largest number of zeros possible)

This one also relates to the linear stuff. It's only valid for particular=True.

>
>
> 'warning=True (default is False)'
> show a warning if checksol() could not conclude.
>
>
> 'simplify=True (default)'
> simplify all but cubic and quartic solutions before
> returning them and (if check is not False) use the
> general simplify function on the solutions and the
> expression obtained when they are substituted into the
> function which should be zero
>
>
> 'manual=True (default is False)'
> do not use the polys/matrix method to solve a system of
> equations, solve them one at a time as you might "manually".
>
> All these parameters trade the speed of the solvers with the accuracy. I
> know there have been issue regarding speed, see
> https://github.com/sympy/sympy/issues/6611. But I still think we should not
> be using these flags for the reason of maintaining simplicity. If a lot of
> people report that our solvers are too slow then we might reintroduce a few
> of them.
>
> It turns out that we might go about not using any of the current parameters
> and have a clean simple input api as
>
> solve(<expression>, <variables>, input_set=<domain of variables>)

I'm not sure if the input_set API is as friendly as it could be. Most
users are going to want to just toss in relations with the rest of the
equations, especially if they can be expressed using the assumptions
or using inequalities, like

solve([sin(x) - 1, x > 0])

or

solve([x**2 + y**2 - z**2, Q.integer(x), Q.integer(y), Q.integer(z)])

Maybe we should just make the API like this. If you want to represent
a complicated set, you can input some kind of Contains object (does
this exist yet?). That would also make it clearer what variables
correspond to what axes in higher dimensional sets.

Aaron Meurer

> https://groups.google.com/d/msgid/sympy/CADN8iurGQmeb5NGD_ZqvB_Y27KQA5h032GmFZRfTGFVEH6n_%3DA%40mail.gmail.com.

[Harsh Gupta]

I'm not sure if the input_set API is as friendly as it could be. Most

users are going to want to just toss in relations with the rest of the
equations, especially if they can be expressed using the assumptions
or using inequalities, like

solve([sin(x) - 1, x > 0])

or

solve([x**2 + y**2 - z**2, Q.integer(x), Q.integer(y), Q.integer(z)])

Maybe we should just make the API like this. If you want to represent
a complicated set, you can input some kind of Contains object (does
this exist yet?). That would also make it clearer what variables
correspond to what axes in higher dimensional sets.

I agree that input_set API is not very friendly but it is consistent for both simple and complicated cases.

I feel that the other methods will get complicated quickly as we try to increase the complexity of the input cases and will add to the maintenance overheads. 

As far as using relations are concerned good conversion tools to set and sufficient examples in documentation shall do as an alternative.

solve(sin(x), x, input_set=(x>0).as_set())

For 

We will have reals as the default parameter for input_set as real so for simple cases user won't have to care about it anyway.

[Harsh Gupta]

Sorry I accidently pushed Send, ( I was unconsciously entering some vim commands). Ignore the "For" after the example.

--
Harsh

[Matthew Rocklin]

> By the way, I didn't know that the sets use exclusively the new assumptions (assuming that is indeed what you are saying here).

That is not what I'm saying.  S.Integers and S.Naturals use new assumptions.  Most of the logic behind core sets doesn't explicitly think about assumptions at all.  So, perhaps worse, sets implicitly uses some combination of the two.  

--

You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CADN8iupYZiJoBNQxCiAtg0Y-bXF1fNs2PbQtdh9qvYd765L7Sw%40mail.gmail.com.

[Matthew Rocklin]

Recall that sets mostly lives in the core.  In particular I think that Interval and maybe Union have been used in the core for some time and FiniteSet has made itself useful in the last couple years.

[Aaron Meurer]

Looking through my notes
(https://github.com/sympy/sympy/wiki/assumptions_handlers, pardon the
formatting; it looks better in emacs), I see that most of the wrong
results I saw were not in the "core" facts, like positive or real
(other than the ones you probably already knew about with infinities
https://github.com/sympy/sympy/issues/5976). It looks like most of the
blatant wrong results that I found were in ntheory (and yes, I am
guilty of not ever opening issues for these things)

In [2]: ask(Q.prime(x*y), Q.integer(x)&Q.integer(y))
Out[2]: False

In [3]: ask(Q.prime(x**y), Q.integer(x)&Q.integer(y))
Out[3]: False

In [6]: ask(Q.integer(abs(x)), ~Q.integer(x))
Out[6]: False

I leave it as an exercise to the reader why each of the above is incorrect.

But more to the point, as I noted in at the bottom of that wiki page,
the new assumptions handlers code is very fragile. If you accidentally
use the wrong kind of Python logic instead of fuzzy logic, you can
return False instead of None. In the best case, you just end up not
implementing as much as you thought you were, but in the worst case,
you get a wrong result.

There are subtle issues that can come up (similar to the
https://github.com/sympy/sympy/issues/5976 one) with the new
assumptions not being consistant in how it defines assumptions. It
doesn't help that it quite often calls out to some external routine to
calculate something, which may or may not use the same definition of
the assumption that it does. It also uses evalf on most assumptions
relating to numbers, which can have its own issues.

Aaron Meurer

> https://groups.google.com/d/msgid/sympy/CAJ8oX-FnZjs4vfBFh%3D2WtMvRdaEM-ndNW%2BqL-rNk5EYQhQ3hVQ%40mail.gmail.com.

[Harsh Gupta]

>> We won't be able to support such output with our new output API
>> assert solve(x - exp(x), x, implicit=True) == [exp(x)]

> How is this implemented?

I read the code and it appears.

Things like exp(x) are replaced by dummy symbols and then solve is allowed to run on it. After the solution is obtained the dummy symbols are replaced back with their original values.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAKgW%3D6L5k4Lvep%3DjVvLPbU6JW%3DMrwsrLsV%2Bq3hAFBQr7XDFdoA%40mail.gmail.com.

For more options, visit https://groups.google.com/d/optout.

--
Harsh

[Harsh Gupta]

We will have to be "very" careful while rewriting the solvers. I just came across another of Professor Fateman's paper which describes why the solution of z**w == y for z isn't simply y**(1/w).

http://www.cs.berkeley.edu/~fateman/papers/y=z2w.pdf

Sadly we doing it wrong.

We are also wrongly simplifying y - (y**(1/(1 + I)))**(1 + I) as 0.

In [58]: f

Out[58]: y - (y**(1/(1 + I)))**(1 + I)

In [59]: simplify(f)

Out[59]: 0

--
Harsh

[Sergey Kirpichev]

On Tuesday, March 18, 2014 8:20:31 PM UTC+4, Harsh Gupta wrote:

We will have to be "very" careful while rewriting the solvers. I just came across another of Professor Fateman's paper which describes why the solution of z**w == y for z isn't simply y**(1/w).

http://www.cs.berkeley.edu/~fateman/papers/y=z2w.pdf

Sadly we doing it wrong.

We are also wrongly simplifying y - (y**(1/(1 + I)))**(1 + I) as 0.

In [58]: f

Out[58]: y - (y**(1/(1 + I)))**(1 + I)

In [59]: simplify(f)

Out[59]: 0

Indeed (e.g. y=-10000 + 4000*I).  Looks a bug somewhere in powsimp().  It doesn't look as a new for me, so maybe it was already reported.

[Aaron Meurer]

If you're working with branch cuts in SymPy, you should be aware of
the very useful exp_polar() function, as well as the related
polar_lift, periodic_argument, and principal_branch.

The basic reason that this is not true is that x**y**z != x**(y*z) for
all x, y, and z, and in particular, for z == 1/y. It is true if x >= 0
or if z is an integer (those are necessary, not sufficient
conditions).

If you haven't studied complex variables, just remember two important things:

- x**y = exp(y*log(x)) (by definition)
- exp(x) == exp(x + 2*pi*I*n) for any integer n

From that (and maybe a little understanding of what log(x) means for
complex x) you can figure out these thing. For instance, x**y**z is
exp(z*log(exp(y*log(x))). If we "cancel" the log and exp, we get
exp(z*y*log(x)), which is x**(y*z). But log(exp(x)) != exp(x) in
general because of the second point, log(exp(x)) == log(exp(x +
2*pi*I*n)). The complex log is multivalued, but we pick the value that
has complex argument in the range (-pi, pi], which gives us a unique
integer n for each complex x such that log(exp(x + 2*pi*n)) == x.

I'll leave it as an exercise to you to combine these and get a
sufficient condition for x**y**z to equal x**(y*z).

You may also be interested in my blog post about a similar problem
http://asmeurersympy.wordpress.com/2013/03/03/when-does-xlogy-ylogx/.

Finally, I should point out that if we are interested in getting *all*
the solutions, then z**w = y -> y**(1/w) is no good anyway. If w is an
integer, there are w solutions (corresponding to the w primitive roots
of unity).

Aaron Meurer

> https://groups.google.com/d/msgid/sympy/CADN8iuocY12bhYgSFkYZCK82jaw1%3D5FGAHFTNoR86ZBuQ9cR9Q%40mail.gmail.com.