I'm willing to invest some of my time to understand if I can be able
to do a step ahead with symbolic functions.
How are special symbolic functions supposed to be defined? I am
willing to experiment with delta of dirac function. This has some
special properties (see http://en.wikipedia.org/wiki/Dirac_delta_function
), some of them are really useful but I don't know how to define them
in a CAS like maxima or SAGE.
I'm aware that it is already present in maxima, even though I don't
think it is recognized by SAGE. I am wondering whether a viable
approach could be to add to calculus.py a section similar to the one
of Function_gamma, so that SAGE simply interfaces to maxima. I don't
know if this is useful or not.
Otherwise, I would be interested in knowing if this could be done with
the new symbolic package. Burcin proved to be very helpful in showing
me a simple way to define delta function by means of its values and he
assigned it as being the derivative of heaviside function (defined in
a sort of piecewise function):
sage: heaviside(x).diff(x)
dirac(x)
is there a way to implement the other properties? I am willing to know
if is there any documentation about that, because I am not able to
find that!
I am willing to learn something about pynac, but please feel free to
discourage me if you think it is too far away from being ready. Is
there any integration or derivation capability ready? Is it possible
to start testing it using maxima's integration capabilities? (I don't
think so...)
I was browsing the todo page ( http://wiki.sagemath.org/symbolics/pynac_todo
) but it seems that many action items went away... are they already
accomplished? (what about the TODO showed in http://wiki.sagemath.org/symbolics
?)
I have to say I find the actual SAGE documentation seems pretty hard
to browse, but could be my fault. Yesterday I lost half an hour
looking for a "numerical solve" or something like that, before finding
the "find_roots" function. Today, I spent half an hour looking for
"differential equation solve" with no success. I am sure I could do
some DE solution in the past (something with maxima, something with
SymPy, I think, all through SAGE), but I think that I found the way to
do it pretty easily browsing the old reference manual...
Thanks a lot
Maurizio
PS: delta of dirac is already in SymPy (
http://code.google.com/p/sympy/issues/detail?id=672&can=1&q=dirac ),
am I correct in thinking that the current function definition is
different than SAGE? In this case, I assume this could be some good
reading, but not necessarily a source of inspiration, right?
On Tue, 14 Apr 2009 14:35:37 -0700 (PDT)
Maurizio <maurizio...@gmail.com> wrote:
>
> Hi all.
>
> I'm willing to invest some of my time to understand if I can be able
> to do a step ahead with symbolic functions.
Great to hear that.
> How are special symbolic functions supposed to be defined?
They will be defined based on the SFunction class in
sage/symbolics/function.pyx. There is no real documentation in that
file, just some examples to test if the wrapper is working as intended.
There is some more documentation in the GiNaC tutorial here:
http://www.ginac.de/tutorial/Symbolic-functions.html
You can do (almost) all that you see in the tutorial, without using any
C++. I say almost, because I still did not set up an equivalent
of .hold(). That should come soon given the current interest in pynac.
Mike Hansen is working on switching Sage's symbolics over to pynac
completely. We talked about changing the design of SFunction on IRC, he
might have some code for you to build on.
Even if you start now, rebasing your code to the new design should be
trivial. You'll just have to put the various functions you define in a
class, which subclasses SFunction.
> I am
> willing to experiment with delta of dirac function. This has some
> special properties (see
> http://en.wikipedia.org/wiki/Dirac_delta_function ), some of them are
> really useful but I don't know how to define them in a CAS like
> maxima or SAGE.
>
> I'm aware that it is already present in maxima, even though I don't
> think it is recognized by SAGE. I am wondering whether a viable
> approach could be to add to calculus.py a section similar to the one
> of Function_gamma, so that SAGE simply interfaces to maxima. I don't
> know if this is useful or not.
With more and more people getting interested in the switch, I think it
has a real chance of happening very soon. I refrained from saying this
before, but at this moment, I don't think it's wise to invest more time
in the old symbolics code.
> Otherwise, I would be interested in knowing if this could be done with
> the new symbolic package. Burcin proved to be very helpful in showing
> me a simple way to define delta function by means of its values and he
> assigned it as being the derivative of heaviside function (defined in
> a sort of piecewise function):
>
> sage: heaviside(x).diff(x)
> dirac(x)
>
> is there a way to implement the other properties? I am willing to know
> if is there any documentation about that, because I am not able to
> find that!
For reference, the code I sent Maurizio with an example containing
heaviside and dirac functions is here:
http://sage.math.washington.edu/home/burcin/pynac/dirac.py
Implementing heaviside is #2452 on trac.
> I am willing to learn something about pynac, but please feel free to
> discourage me if you think it is too far away from being ready. Is
> there any integration or derivation capability ready? Is it possible
> to start testing it using maxima's integration capabilities? (I don't
> think so...)
Derivation is easy, and it's already done by pynac. We will be using
maxima for integration for a while now, though I plan to start
submitting native integration (and summation) code to Sage once the
switch is complete.
For people interested in helping out, a few words about integration is
in order. As we discussed on this list before, many problems in
integration can be handled with heuristics before calling more advanced
algorithms. The heuristics also give more user friendly answers in most
cases.
Integration code in Sage should run through some preprocessing before
calling anything else, be it maxima, or a native implementation of the
Risch algorithm.
There are some links to descriptions of heuristics used by maxima and maple at this page:
http://wiki.sagemath.org/symbolics
Pattern matching and rewriting capabilities of the new symbolics will
be very useful for this. Anybody interested in improving the
integration capabilities of Sage can help, no high level mathematical
knowledge is necessary.
> I was browsing the todo page
> ( http://wiki.sagemath.org/symbolics/pynac_todo ) but it seems that
> many action items went away... are they already accomplished?
The items in the todo list didn't go away, they were copied to the
"done" section, with a link to the trac ticket containing the changes.
There are two pynac packages waiting for review on trac. Since they
both contain substantial changes, it would be good if somebody other
than the developer took them for a ride. :)
I only know two people who tried to use pynac for anything real until
now, Alex Raichev and Nick Alexander. Their comments have been really
helpful and motivating.
> (what
> about the TODO showed in http://wiki.sagemath.org/symbolics ?)
This is a more long term todo list, both to remember what still needs
to be done and for reference if anybody is interested in helping out.
Most tasks on this list are highly nontrivial, and need substantial
work to be complete.
Feel free to edit the list and add items you think are important. For
example, you can add z-transforms and laplace transforms as we
discussed before.
Problems related to special functions used to be classified under the
component calculus on trac. I am not sure if I should add special
functions stuff to that page, or start a new one. We have quite a few
special functions related bugs on trac. There was also some effort to
get Sage on the list of software in the DLMF, but the wiki page for
this seems to be deleted. (mabshoff?)
> I have to say I find the actual SAGE documentation seems pretty hard
> to browse, but could be my fault. Yesterday I lost half an hour
> looking for a "numerical solve" or something like that, before finding
> the "find_roots" function. Today, I spent half an hour looking for
> "differential equation solve" with no success. I am sure I could do
> some DE solution in the past (something with maxima, something with
> SymPy, I think, all through SAGE), but I think that I found the way to
> do it pretty easily browsing the old reference manual...
sage: deso<tab>
desolve desolve_laplace desolve_system desolvers
Especially "desolvers?" seems helpful.
> Thanks a lot
Thank you.
Cheers,
Burcin
No problem. We didn't lose anything much, since the "effort" didn't get
anywhere. The page Robert Bradshaw started contained a simple header,
and a few section headings without any real content. It won't be hard
to reproduce it, if there is any time that is...
BTW, could we get the subscriptions for the wiki working again?
Thanks.
Burcin
A is 1 in example 1.5.
> B = x^3 + x
> tores = A - z*B.derivative(x)
> res = tores.resultant(B,x); factor(res)
>
> I get:
> (-1) * (z + 1) * (2*z - 1)^2
With A = 1:
sage: P.<x,z> = QQ[]
sage: b = x^3+x
sage: b.resultant(1-z*b.derivative(x),x)
-4*z^3 + 3*z + 1
sage: factor(_)
(-1) * (z - 1) * (2*z + 1)^2
which matches the maple result.
> That is somehow close to the result shown in the link, but why is it
> not identical?
>
> Moreover, why if I use the second statement (the one in the comment)
> to define the ring, do I get a completely different answer? E.g.:
> (z + 1) * (z + 2)^2
>
> Which is the correct way to define the ring in this case?
The theory only works over characteristic 0, i.e., your fields should contain QQ. Also note that,
sage: P.<x,z> = GF(5)[]
sage: (x+x^5).derivative(x)
1
What is the integral of 1 w.r.t x in this case?
> Finally, even assuming that I can get the right answer from this,
> which is the recommended way to get the roots of an equation given by
> a "univariate polynomials == 0"? This is supposed to be the next step
> of the algorithm.
>
> The worst way I can think of is to use repr(), so that I can use solve
> (), but hopefully there's a better way, at least to convert a
> polynomials to a symbolic expression... Is there?
You can get this information from the factorization you compute above.
Cheers,
Burcin
Taking a quick look at that page, it looks like they want the exact
roots in CC of a polynomial with algebraic coefficients. In Sage, we
can get this with QQbar:
sage: K.<x> = QQbar[]
sage: (x^5-x-1).roots(ring=QQbar)
[(1.167303978261419?, 1),
(-0.7648844336005847? - 0.3524715460317263?*I, 1),
(-0.7648844336005847? + 0.3524715460317263?*I, 1),
(0.1812324444698754? - 1.083954101317711?*I, 1),
(0.1812324444698754? + 1.083954101317711?*I, 1)]
(AFAIK, maxima can't do this; I don't think maxima can handle general
algebraic numbers. Since Sage's solve() is implemented using maxima,
solve() won't work for this problem.)
Carl
GF(5) means the Galois field of characteristic 5, a.k.a. the integers
modulo 5. (So in GF(5), you have 2*3 = 1, 3+4 = 2, etc.) It's
probably quite irrelevant for computing integrals.
QQ is the rational numbers (fractions). QQbar is the algebraic
closure of QQ; this means it includes every complex number which is
the root of a polynomial with rational coefficients. So it includes
things like sqrt(2) (which is a root of x^2-2), and sqrt(-1) (a root
of x^2+1), as well as more exotic numbers like the roots of x^5-x-1,
which can't be expressed using radicals (roots). (QQbar does not
include all complex numbers, though; for instance, it does not include
pi or e, which are transcendental rather than algebraic.)
> Now let's go to Carl's help...
>
>> Taking a quick look at that page, it looks like they want the exact
>> roots in CC of a polynomial with algebraic coefficients. In Sage, we
>> can get this with QQbar:
>>
>> sage: K.<x> = QQbar[]
>> sage: (x^5-x-1).roots(ring=QQbar)
>>
>
> First problem with QQbar: it seems that resultant() doesn't like it,
> because it is not able to convert it to a Singular ring (this is the
> error, I'm not attaching all the output, tell me if you need it)
>
> TypeError: no conversion of this ring to a Singular ring defined
Looks like this hasn't been implemented yet.
> On the contrary, QQ[] seems to work fine with resultant (but it
> doesn't have roots() )
Univariate polynomials over QQ definitely have roots(); were you using
a multivariate polynomial ring?
> Moreover, it seems that QQbar roots() is not working for multivariate
> polynomials ring... is it true or am I just missing something else? In
> that case, is possible to let it work in multivariate polynomials? As
> you can imagine, I would like to think about this as a method of
> solving integrals, so it is very likely to have a symbolic expression
> with more than just a single symbolic variable.
>
> Apart from this, is there another way to solve an equation (with more
> than a single symbolic variable) obtaining exact roots? It seems that
> maxima would do the work (with algebraic numbers...), is it possible
> that it is the only symbolic equation solver within SAGE? What about
> SymPy or anything else?
What does this even mean? .roots() gives a list of all the solutions
of a univariate polynomial equation. But a multivariate polynomial
equation will usually have an infinite number of solutions; for
instance, x^2+y^2-1=0 has an infinite number of solutions over the
rationals (or the reals, or the algebraic reals, etc.)
If you have a system of multivariate polynomial equations, then the
system might have only finitely many solutions.
> Finally, I still would like to know which is the best way to translate
> the output of a calculation with polynomial rings into a symbolic
> expression, that can be carried on with maxima or pynac. Can you help
> me?
If p is a polynomial, then SR(p) is a symbolic expression.
Carl
It's not elegant, but it's probably much more efficient... QQbar is
vastly less efficient than QQ, so it's definitely a good idea to stick
with QQ as long as possible before jumping to QQbar. You can use
.roots(ring=QQbar) on a polynomial over QQ.
> if I don't define a ring containing x and z, how do I inject z into
> the workspace? certainly not as a symbolic variable, so what else? I
> know that at the very end this is a univariate polynomial ring root
> finding problem, I'm just wondering how to correctly manage it.
If you have a polynomial in a multivariate ring, but your actual
polynomial actually contains only one variable, you can get the
corresponding univariate polynomial with .univariate_polynomial().
> I understand, but let's consider the problem I'm coming from:
> indefinite integral. I'm not using polynomials because I'm well aware
> of the implications, rather because I found resultant() implemented
> there :)
> What I intended with my sentences was: imagine a polynomial like: x^2
> + a^2. The solutions of x^2 - a^2 = 0 with respect to x over any field
> wouldn't always coincide with +/- a? That is the kind of reasoning I
> would make to solve the indefinite integral, like when I work with
> symbolic variables that have a precise meaning (i.e. they are not
> really variables in the problem, it's just that I want to keep that
> degree of freedom with respect to the solution, especially if I want
> to see the effect of that value over the solution) into my problem.
I think it's going to be almost impossible to carry through the
algorithm with parameters... it's definitely not where you should be
starting :)
Consider x^5-x-a. In general I don't know of a better description of
the roots of this polynomial as a function of a than "the roots of
x^5-x-a as a function of a" :) (And if you want the distinct roots,
that's even worse. The number of distinct roots will depend on the
values of the parameters. Fortunately it probably doesn't matter
whether the roots are distinct, if you carefully adapt the
algorithm...)
>> > Finally, I still would like to know which is the best way to translate
>> > the output of a calculation with polynomial rings into a symbolic
>> > expression, that can be carried on with maxima or pynac. Can you help
>> > me?
>>
>> If p is a polynomial, then SR(p) is a symbolic expression.
>>
>
> Funny, the next minutes I spent, I recognized how to do that in this
> same way :) By the way, is there a SR equivalent for the new pynac
> symbolic? Something like NSR (new symbolic ring)? I can see that pynac
> already has gcd implemented, although I can already (unfortunately)
> see:
> RuntimeError: gcd: arguments must be polynomials over the rationals
There is an NSR, but it's not available on the command line by
default; you can get it with:
from sage.symbolic.ring import NSR
Carl
Well, we just need a resultant algorithm that doesn't go through
Singular. I'm planning to write such a thing as part of my
cylindrical algebraic decomposition implementation sometime in the
next few months.
Carl
Looks like a bug to me.
Burcin, any comments?
> Moreover, I don't see this being the right way to do this, because
> (for this particular problem: integration) I don't like having the
> numerical representation of things like sqrt(5), even if the result is
> still correct, so that
>
> temp = QQbar(sqrt(5)); temp
> 2.236067977499790?
>
> temp^2
> 5.000000000000000?
Note that it's only the representation that's numerical; internally
these are actually exact numbers. (For instance, if you do "temp^2 ==
5" and Sage prints "True", then that means that temp really is exactly
the square root of 5.) I use a numeric representation because not all
algebraic numbers can be represented with radicals, and solving by
radicals is much more difficult than the symbolic-numeric algorithms I
used in QQbar.
> So, please tell me. Which should be the right way to try to approach
> this indefinite integration problem? You can see that I'm not that
> good in deep mathematical theory, but approaching the simplest problem
> (that could be different from this one I'm looking at right now) is
> fun :)
Unfortunately, I think the right way is a lot of work.
1) Teach yourself enough math that you can mostly understand the
algorithm in the paper. (Not just mechanically follow the steps, but
understand what each step does, and have some idea why the algorithm
as a whole works.)
For some parts of this task, wikipedia and mathworld.wolfram.com will
be useful resources (planetmath.org is also interersting); but you'll
probably also need to read some books, etc.
2) Implement each missing sub-algorithm in Sage.
3) Implement the entire algorithm in Sage.
This sounds glib, but it really is possible... I'm implementing
cylindrical algebraic decomposition in Sage following this program.
I've been working on it in my spare time for several years; I probably
spent a couple of years on step 1, then I started step 2 a couple of
years ago, and now I'm working on step 2 and step 3.
The problem is that the Risch algorithm just is not simple.
A good integration algorithm probably has a couple of phases -- first
you try some heuristics, and pattern-match against a database of known
integrals; if those fail, then you bring out the big guns and apply
the decision procedure. If you're interested in working on
integration, maybe you'd rather work on the first phase? That
probably requires more computer science and less math than the second
phase.
Carl
Curiously though, SymPy knows this particular integral.
>>> from sympy import *
>>> x=Symbol('x')
>>> f=(x**2+2*x+1+(3*x+1)*sqrt(x+log(x)))/(x*sqrt(x+log(x))*(x+sqrt(x+log(x))))
>>> integrate(f,x)
2*log(x + (x + log(x))**(1/2)) + 2*(x + log(x))**(1/2)
Fredrik
A much shorter example is:
integrate(sqrt(x+log(x)),x)
to which Axiom replies:
integrate: implementation incomplete (constant residues)
Tim Daly
A much shorter example is:
integrate(sqrt(x+log(x)),x)
to which Axiom replies:
integrate: implementation incomplete (constant residues)
For further information, Bronstein deals with constant residues
(written after his Axiom work) in his book [1] after page 147.
Tim Daly
[1] Bronstein, Manuel "Symbolic Integration I", Second Edition
Springer ISBN 3-540-21493-3
What is f(x) = sqrt(x+log(x)) supposed to be an example of? Does f
has an antiderivative that can be expressed in terms of elementary
functions? If so, what is it?
-- William
I was just about to point it out. It even works in Sage:
In [2]: f = (x^2 + 2*x + 1 +
...: (3*x+1)*sqrt(x+log(x)))/(x*sqrt(x+log(x))*(x+sqrt(x+log(x))))
In [5]: import sympy
In [6]: sympy.integrate(f, x)
Out[6]: 2*log(x + (x + log(x))**(1/2)) + 2*(x + log(x))**(1/2)
So we have a good start to implement the Risch algorithm in sympy already.
Ondrej
Ondrej,
How does sympy compute the integral of
(x^2 + 2*x + 1 + (3*x+1)*sqrt(x+log(x)))/(x*sqrt(x+log(x))*(x+sqrt(x+log(x))))?
William
Ondrej, what result do you get for:
integrate(sqrt(x+log(x)),x)
Tim
Axiom has one of the most complete Risch implementations. Many of
the original authors of the theory worked on the implementation.
Axiom can almost do the Risch integral of the formula you posted.
It fails with:
integrate: implementation incomplete (constant residues)
The failing subexpression sqrt(x+log(x)) falls into the unimplemented
case. Replacing this subexpression by 'a' succeeds.
Most systems do pattern matching and heuristics. Fateman had a server
available that would "integrate" using these techniques.
Tim
In [1]: integrate(sqrt(x+log(x)),x)
Out[1]:
⌠
⎮ ⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽⎽
⎮ ╲╱ x + log(x) dx
⌡
So we can't do it yet.
Ondrej
> Ondrej,
>
> How does sympy compute the integral of
>
> (x^2 + 2*x + 1 + (3*x+1)*sqrt(x+log(x)))/(x*sqrt(x+log(x))*(x+sqrt(x+log(x))))?
See this thread in our list:
http://groups.google.com/group/sympy/browse_thread/thread/9a61d681e6f967e8/
Ondrej
>
> On Sun, Apr 19, 2009 at 7:44 AM, Maurizio
> <maurizio...@gmail.com> wrote:
> > Carl, I took advantage of your suggestion, even though I assume I
> > can't still go through the whole process with the current gcd
> > capabilities in Pynac. But before than that, I'd like to point out
> > something strange I did notice, and maybe also Burcin can help with
> > that:
> >
> > reset()
> > P.<x,z> = QQ[]
> >
> > B = x^3 + x
> >
> > var('x, zs', ns = 1)
> > from sage.symbolic.ring import NSR
> > Bs = NSR(B)
> > Bs
> > x^3 + x
> > Bs.diff(x)
> > 0
> >
> > So, the derivative is not working. Which is the cause? It seems that
> > the "x" in Bs is not the "x" I declared, so the derivative gets 0
> > as a result. Which is the reason?
>
> Looks like a bug to me.
>
> Burcin, any comments?
I agree that it's confusing, but it's not a bug.
The command
sage: Bs = NSR(B)
converts the polynomial B = x^3 + x in QQ[x] to a symbolic expression,
with one numeric coefficient, namely B.
sage: Bs.coeff(x)
0
sage: Bs.is_constant()
True
sage: Bs.pyobject()
x^3 + x
sage: Bs.pyobject().parent()
Multivariate Polynomial Ring in x, z over Rational Field
Since Bs is a constant w.r.t. the symbolic x, the derivative of Bs
w.r.t. the symbolic variable x is 0.
In order to convert the polynomial to a symbolic expression, I suggest
substituting the symbolic variables in the polynomial.
sage: B.subs(x=x)
x^3 + x
sage: t = B.subs(x=x)
sage: t.coeff(x)
1
sage: t.coeff(x^3)
1
sage: t.derivative(x)
3*x^2 + 1
<snip>
> A good integration algorithm probably has a couple of phases -- first
> you try some heuristics, and pattern-match against a database of known
> integrals; if those fail, then you bring out the big guns and apply
> the decision procedure. If you're interested in working on
> integration, maybe you'd rather work on the first phase? That
> probably requires more computer science and less math than the second
> phase.
Indeed, this was what I tried to say in those comments about
integration that got Maurizio started on this track.
There are descriptions of heuristics we can use in the maxima manual [1]
and the book by Geddes, Czapor and Labahn [2].
[1]
http://maxima.sourceforge.net/docs/manual/en/maxima_20.html#Item_003a-integrate
[2] http://books.google.at/books?id=9fOUwkkRxT4C&pg=PR2&hl=en#PPA473,M1
It would be great if someone went through these links, and made them
more explicit/precise. As I said before, implementing these should be
within reach with the pattern matching/rewriting capabilities of pynac.
Here is an example, pointed out to me by Clemens Raab, which maxima can
do and MMA can't:
sage: g=log(1-x)/(polylog(2,x)^2+log(1-x)*polylog(3,x))
sage: f = g.derivative(x)
sage: integral(f,x)
log(1 - x)/(log(1 - x)*polylog(3, x) + polylog(2, x)^2)
Cheers,
Burcin
So is this what happens anytime you do NSR(sage_object)? It converts it
to a constant (as far as symbolics is concerned) with a coefficient of
sage_object? Or are there any exceptions?
Jason
--
Jason Grout
I consider this bad behavior as well. Perhaps it should at least try
to call a _symbolic_ method.
- Robert
> Hi Michael,
>
> Actually, I thought that this discussion (especially people much more
> expert than me) has clarified the point that implementing integrals is
> not really just matter of a couple of months... but I would be glad to
> see this happen!
Implementing something that can handle a huge number of integrals is
a reasonable goal for a couple of months. Implementing something that
can handle everything that we know how to handle in theory...well,
that hasn't ever happened yet.
> I know there are some license issues with SymPy (not really issues,
> just differences probably) but I think there's no problem in taking
> inspiration and some pieces of code from it, right?
There is a problem with just lifting code--but we can and do ship
SymPy as part of Sage.
> I'm saying this, because I can see this new symbolic stuff exciting,
> but without the right amount of interest on it, its development will
> always be a little slow
Given that we ship SymPy, so anything it can handle we should be able
to handle. I imagine when you want to integrate something, it will
try to do it natively first, and that failing ask SymPy and/or maxima
before returning a failure. Over time we'll depend less and less on
the other two.
- Robert
They are --- you can use them from sympy inside Sage. It's my goal
that all sympy features are nicely integrated in Sage. I work on this
as time permits.
> At least, they are already aware of their shortcomings (ie: cannot
> compute the integral of log(x)/x ).
>
> I'm sure SAGE people could give big contribute to those, send patches
> upstream, and move forward if needed.
We got 5 summer of code students this year, so we'll move sympy a lot
forward this summer. The most important project:
http://socghop.appspot.com/student_project/show/google/gsoc2009/portland_state/t124024247737
will disentangle the assumptions, so that we should then be able to
easily use any other core, like pynac, or our own cython core. Besides
that when my school ends in less than a month, I plan to work hard on
SPD:
http://code.google.com/p/spdproject/
so that one can easily create Sage like projects with his own
packages. I expect this to make it easier for people to create
libraries in a way, so that they are easy to use with Sage, e.g. they
work with plotting, etc., but they could distribute their own version
of all in one distribution.
In any case, this will be an exciting summer.
Ondrej
Also, in the pynac-based symbolics that Mike Hansen has been polishing
up for full inclusion in Sage (to replace the maxima based symbolics),
one can just do
f.integrate(algorithm="sympy")
and sage will compute the integral using sympy. He's already
implemented this and I've seen it work well when I tried it out.
>> At least, they are already aware of their shortcomings (ie: cannot
>> compute the integral of log(x)/x ).
>>
>> I'm sure SAGE people could give big contribute to those, send patches
>> upstream, and move forward if needed.
>
> We got 5 summer of code students this year, so we'll move sympy a lot
> forward this summer. The most important project:
>
> http://socghop.appspot.com/student_project/show/google/gsoc2009/portland_state/t124024247737
>
> will disentangle the assumptions, so that we should then be able to
> easily use any other core, like pynac, or our own cython core. Besides
> that when my school ends in less than a month, I plan to work hard on
> SPD:
>
> http://code.google.com/p/spdproject/
>
> so that one can easily create Sage like projects with his own
> packages. I expect this to make it easier for people to create
> libraries in a way, so that they are easy to use with Sage, e.g. they
> work with plotting, etc., but they could distribute their own version
> of all in one distribution.
>
> In any case, this will be an exciting summer.
Cool. Are you going to keep Microsoft Windows in mind wrt SPD? I
think the work Dan Shumow (and others) have been
doing on the Sage windows port could, in a sense, be seen as a base
for a completely open SPD for Windows. See
http://windows.sagemath.org/
-- William
Yes, that way Sage allows to call any of those libraries very easily
in Sage and at the same time uses just one library by default, that
currently is the best at integrating, I guess still Maxima.
>
>>> At least, they are already aware of their shortcomings (ie: cannot
>>> compute the integral of log(x)/x ).
>>>
>>> I'm sure SAGE people could give big contribute to those, send patches
>>> upstream, and move forward if needed.
>>
>> We got 5 summer of code students this year, so we'll move sympy a lot
>> forward this summer. The most important project:
>>
>> http://socghop.appspot.com/student_project/show/google/gsoc2009/portland_state/t124024247737
>>
>> will disentangle the assumptions, so that we should then be able to
>> easily use any other core, like pynac, or our own cython core. Besides
>> that when my school ends in less than a month, I plan to work hard on
>> SPD:
>>
>> http://code.google.com/p/spdproject/
>>
>> so that one can easily create Sage like projects with his own
>> packages. I expect this to make it easier for people to create
>> libraries in a way, so that they are easy to use with Sage, e.g. they
>> work with plotting, etc., but they could distribute their own version
>> of all in one distribution.
>>
>> In any case, this will be an exciting summer.
>
> Cool. Are you going to keep Microsoft Windows in mind wrt SPD? I
Yes, definitely, working on windows is one of the goals.
> think the work Dan Shumow (and others) have been
> doing on the Sage windows port could, in a sense, be seen as a base
> for a completely open SPD for Windows. See
> http://windows.sagemath.org/
Excellent, I just joined the list and will work closely with him.
Ondrej
unless maxima doesn't produce an answer (i.e., returns your expression
unevaluated or throws an exception because of an assumption, etc.) Then
it could automatically try sympy.
-Jason
--
Jason Grout
Excellent. Please report any bugs that you discover in the sympy<->
sage conversion. I believe there still may be some.
Thanks,
Ondrej
Please don't report any bugs you discover in the sympy<-->sage conversion
right now, since Mike Hansen literally wrote a bunch of it a few days ago,
and it exists only on his laptop.
Once Sage-4.0 appears with the new symbolics, then we'll definitely
want to hear about bugs.
William
--
William Stein
Associate Professor of Mathematics
University of Washington
http://wstein.org
I believe that too.
>
>> > I know there are some license issues with SymPy (not really issues,
>> > just differences probably) but I think there's no problem in taking
>> > inspiration and some pieces of code from it, right?
>>
>> There is a problem with just lifting code--but we can and do ship
>> SymPy as part of Sage.
>
> Yes. One could take code from Sympy (the BSD license allows this) and
> make GPL ed changes on top of it. The main issue is just that our
> pattern matching engine via pynac will be more powerful, our
Could you give some specific examples where the pynac pattern engine
is more powerful?
> arithmetic is faster and the other abstract math bits in Sage are way
> more powerful than what is in Sympy. And the goal of Sympy is to be
> self contained and depending on certain operations in Sage is not an
> option and not wanted, so contributing such code back to Sympy is not
We want sympy to run without Sage. But I have absolutely no problems
calling Sage for certain things if it's available. In fact, I do want
to call Sage if it's available.
> an option . There is the goal for Sympy to optionally depend on Pynac
> and use it when available, but no one is working on this, so this is
> not something we ought to be waiting for.
>
> Plainly and simply: Waiting on someone else to fix the problem for us
> has not worked, so we will do it ourselves.
For the record, I offered William and you that I will work on this, if
we manage to find funding for it over the summer. So just to be clear
that I am very interested in this, and not just talking. But
unfortunately I myself didn't find funding for it (and I haven't heard
from you either that you found some possibility), so I had to find an
internship somewhere else related more to my research (finite
elements).
We managed to get one gsoc project that does the assumptions right, so
it may happen anyways over the summer, in fact I very much hope so.
>> > I'm saying this, because I can see this new symbolic stuff exciting,
>> > but without the right amount of interest on it, its development will
>> > always be a little slow
>>
>> Given that we ship SymPy, so anything it can handle we should be able
>> to handle. I imagine when you want to integrate something, it will
>> try to do it natively first, and that failing ask SymPy and/or maxima
>> before returning a failure. Over time we'll depend less and less on
>> the other two.
>
> Yep, that is the only way in my opinion. We control the Sage library
> and can coordinate releases with improvements in the Sage library -
> Sympy has not released as frequently as it used to be and unless
> someone steps up and helps Ondrej out I doubt that is going to change
> any time soon.
We'll see. We managed to get 5 gsoc student this summer, so we are not
at all dead, if this is what you mean. :)
Besides, it will take almost a year to Sage too to release the new
symbolics (started August 2008), so I don't think we are doing that
bad.
Also in terms of developers working just on symbolics, I think sympy
has much more developers. I don't know if it's easy to extract just
patches to Sage symbolics, to compare speeds of developments, but I
would not be surprised if it turns out it's actually very comparable,
or even less people are working on Sage symbolics than on sympy.
That being said, I like that you pursue the pynac way, because first I
also wanted to use ginac but using swig was just not the way, so
William showed me with cython how to actually do it. And second, I
welcome competition, because that's the only way to actually move
forward, but for Sage and sympy. For example thanks to sympy, you
completely abandoned your previous approach for symbolics, because
after many months of development, it was still slower than pure python
sympy and less tested and robust. So now it's our turn to speed up
sympy. :)
Ondrej
but -> both
O.
There is an integration test suite at:
http://axiom-developer.org/axiom-website/CATS
which has the Schaums integral series along with examples.
Each integral result is subtracted from the Schaums answer
and then simplified, hopefully to a constant.
How does Sympy compare?
Tim Daly
As far as I know, it hasn't been tested with it yet. I'm planning support for
it, but I need to do some work in order to get the test suite working with it.
It's on my list.
Cheers,
Tim Lahey.
It's on our list too, so it will happen eventually. We definitely
still need to improve our algorithms a lot, see e.g.:
http://groups.google.com/group/sympy/browse_thread/thread/58916fb31e1ff1ea
but a nice thing is that it's in Python, so it's easy to work with.
Ondrej
At some point I probably should get around to joining that group.
If you have any suggestions for the testing architecture, I'd be happy to hear
them.
I've now started using GitHub's issues to do bug/issue tracking for the code.
Cheers,
Tim Lahey.
We use py.test/nosetest compatible tests, but if you prefer Sage like
doctests, that's fine too.
> I've now started using GitHub's issues to do bug/issue tracking for the code.
I only use github for storing the git repository, haven't used their
bugs/issues yet.
Ondrej
PDE's are of course important, also it's my main research thing, but
for sympy the assumptions are essential, because you cannot build a
more advanced symbolics without a good assumption system. I am
curious which approach Sage is going to use for this, since ginac
doesn't have any assumptions.
>>
>> We'll see. We managed to get 5 gsoc student this summer, so we are not
>> at all dead, if this is what you mean. :)
>>
>
> That's definitely a very good thing, especially in getting people
> involved with development... Probably the best result is not just
> their short term, but also the chance of some long term commitment of
> good people.
Exactly. Getting people to work with you is absolutely essential,
that's my priority number one. If you work alone, you cannot do
anything in the long term.
Ondrej
Good to know. I'm not actually using any unit test classes at the moment
since I have the difficulty that each test has different steps to properly
compare the output to the Schaum's result. The multiple solutions also
pose an additional difficulty.
Cheers,
Tim Lahey.
I think sympy will do very poorly if the assumptions are needed, we
are still working on the assumptions.
So if it turns out too difficult, just skip sympy for the time being,
we'll get back to it later.
Ondrej
The problem arises with all the different integration systems. Usually some
kind of simplification is needed on the integral returned, even if there aren't
multiple solutions. This complicates the testing procedure since the steps to
perform the simplification are often specific to the returned result.
I'm currently aiming to finish the test suite just for Sage/Maxima and
I'll go back
and address the various issues (like testing SymPy) once that's complete.
Cheers,
Tim Lahey.
Ok, awesome. I am interested in all bugs that you find.
Thanks for all the work you are doing,
Ondrej
Would it be better to test the results numerically? (For instance,
evaluate the integral returned and the desired result at 100 random
points to high precision, and ensure that the relative error between
the answers at each point is small.)
Of course, this wouldn't count as a proof that the result was correct,
but IMHO it would be good enough (it seems unlikely that integration
bugs would result in wrong answers that were numerically almost
equivalent to the right answer). (Actually, I might actually trust a
numeric result more than a symbolic simplification-based result, given
the theoretical possibility that a simplification bug might cancel out
an integration bug, leading to a false pass in the test suite;
especially if simplification and integration are done in the same
system.)
Carl
How can one do that with symbolic variables? Most of the test integrals
have symbolic constants (w.r.t the integration) so it isn't just the
integration
variable. I thought about numerical testing, but it isn't generally feasible.
>
> Of course, this wouldn't count as a proof that the result was correct,
> but IMHO it would be good enough (it seems unlikely that integration
> bugs would result in wrong answers that were numerically almost
> equivalent to the right answer). (Actually, I might actually trust a
> numeric result more than a symbolic simplification-based result, given
> the theoretical possibility that a simplification bug might cancel out
> an integration bug, leading to a false pass in the test suite;
> especially if simplification and integration are done in the same
> system.)
It's possible to have simplification bugs, but I'll have to rely upon separate
tests of the simplification system.
Cheers,
Tim Lahey.
Couldn't you just pick random values for all of the symbolic constants, as well?
Carl
Yes, but over what range? If you do that, you've just ensured that it
is correct for
those points. It also could get expensive if you have multiple
constants combined
with the integration variable. It's an option to consider, but I'd
probably only consider
for those that don't simplify to the same result. It could also be
done as an optional
flag. I'll add it to the list of things to consider.
If someone wants to add issues/bugs for me, it can be done at:
http://github.com/tjl/sage_int_testing/issues
It requires a GitHub account (which is free). Or one can e-mail me
directly and I'll
add it.
Cheers,
Tim.
What's the github name?
Tim Daly
The issue tracker:
http://github.com/tjl/sage_int_testing/issues
and the main repository:
http://github.com/tjl/sage_int_testing/tree/master
Cheers,
Tim.
---
Tim Lahey
PhD Candidate, Systems Design Engineering
University of Waterloo
http://www.linkedin.com/in/timlahey
Is there anyone working on an assumptions framework in Sage? On the one
hand, you can work in certain domains in Sage (i.e., polynomials over
QQ, etc.), so some of this need may be taken care of there.
But for a more general framework (like declaring that x>0?)
Hmmm...maybe we'll use sympy? :)
Jason
--
Jason Grout
>
> Could you explain how assumptions are so important? Could you
> particularly address how they can (1) be so critically important, and
> yet (2) ginac doesn't have them. Incidentaly, to me they are
> particularly important in symbolic integration, which ginac doesn't
> do. Also, could you explain why the assumption system in Sympy
> needs to be rewritten -- in particular, what design decisions were
> suboptimal?
>
> -- William
>
One place where I use assumptions regularly is with trig functions. If
you have sin(n*pi) or cos(n*pi), stating the assumption that n is
integer
reduces the first to 0 and the second to alternating +1,-1. These crop
up
in modal analysis of physical systems regularly. I don't know
if sympy's assumptions are suboptimal, since I haven't used sympy much,
but I guess it is since they're planning on improvements.
We already discussed this many times on this list, just search the
archives. Without good assumptions, you cannot implement good
integration, good solvers, good simplification, nothing. Of course,
things like x+x is easy, but once you have for example
Integral(<complex expression involving many constants>), sometimes you
can simplify it, sometimes not and this very much depends on the
assumptions for the constants. Or things like integrate(x**n, x).
>> particularly address how they can (1) be so critically important, and
>> yet (2) ginac doesn't have them. Incidentaly, to me they are
ginac doesn't have any assumptions and so ginac doesn't have any
advanced symbolic features.
Of course, if the only thing that you are going to use ginac for is
the core engine, then it should not matter much. But what I want from
a CAS is to be able to do advanced symbolic calculations with symbols,
so I need some way to tell it my assumptions about the symbols. Taking
everything as complex numbers will not simplify things in many cases.
>> particularly important in symbolic integration, which ginac doesn't
>> do. Also, could you explain why the assumption system in Sympy
>> needs to be rewritten -- in particular, what design decisions were
>> suboptimal?
Because we attach assumptions to symbols, e.g. you define
In [2]: x = Symbol("x", positive=True)
and then you work with that:
In [3]: sqrt(x**2)
Out[3]: x
In [4]: x = Symbol("x", negative=True)
In [5]: sqrt(x**2)
Out[5]: -x
In [6]: x = Symbol("x")
In [7]: sqrt(x**2)
Out[7]:
⎽⎽⎽⎽
╱ 2
╲╱ x
But this approach is the wrong one, because then the core engine has
to take this into account and it slows things down. Our new system
will keep the assumptions external, and define methods to work with
them, like in mathematica, e.g. refine(). This should simplify the
core and then we should be able to speed it up.
>>
>> -- William
>>
>
> One place where I use assumptions regularly is with trig functions. If
> you have sin(n*pi) or cos(n*pi), stating the assumption that n is
> integer
> reduces the first to 0 and the second to alternating +1,-1. These crop
> up
> in modal analysis of physical systems regularly. I don't know
> if sympy's assumptions are suboptimal, since I haven't used sympy much,
> but I guess it is since they're planning on improvements.
Exactly, all kinds of these assumptions, that are necessary for
integration, simplifications, e.g. sometimes sqrt(x**2) reduces to x,
sometimes to -x, sometimes to neither.
Ondrej
Note that ginac has a similar mechanism, which it uses for the
automatic simplifications it does. I haven't exposed this in the Sage
interface. I don't know if Mike has either.
Here is a list of the info flags ginac supports (Scroll down for the
table):
http://www.ginac.de/tutorial/Information-about-expressions.html
Pynac has an "infinity" flag in addition to the ones listed there.
I admit that I have no experience with general simplification
routines, but I have the feeling that it all boils down to answering
is_positive(), is_integer() queries about expressions. The trick is to
be able to deduce the answers to these from a few bits of information
gathered from the user.
I don't know of any open source package which can provide this
capability for Sage. I've heard that Reduce with its Redlog package is
very powerful in this area. I doubt if it can be used from Sage though.
Cheers,
Burcin