[PATCH] [mq]: integration.patch
Ondrej Certik
> integration.patch
> 9KDownload
Actually, I made a mistake here:
diff -r 315e41f3d28d -r 42b943e63382 examples/fem.py
--- a/examples/fem.py Tue Nov 27 01:38:50 2007 +0100
+++ b/examples/fem.py Tue Nov 27 02:28:12 2007 +0100
@@ -30,7 +30,7 @@ class ReferenceSimplex:
for d in range(0,nsd):
p = coords[d]
limit += p
- intf = integrate(intf, (p, 0, limit))
+ intf = integrate(intf.expand(), (p, 0, limit))
return intf
def bernstein_space(order, nsd):
The above change shouldn't be present in the patch. Fortunately we use
a patch review now. :)
I delete the .expand() in the example and it still executes in 2.4s,
because it is expanded while converting to a polynomial. So my patch
takes the original example and speeds it up by a factor of 10. Not
bad.
Kirill Smelkov
> # HG changeset patch
> # User Ondrej Certik <ond...@certik.cz>
> # Date 1196126892 -3600
> # Node ID 42b943e63382333c66eafdaf2d396bb5b07b7aa2
> # Parent 315e41f3d28d0c7bd505b86da8a58d9ebb48288e
> [mq]: integration.patch
Log message would be nice.
>
> diff -r 315e41f3d28d -r 42b943e63382 examples/fem.py
> --- a/examples/fem.py Tue Nov 27 01:38:50 2007 +0100
> +++ b/examples/fem.py Tue Nov 27 02:28:12 2007 +0100
> @@ -30,7 +30,7 @@ class ReferenceSimplex:
> for d in range(0,nsd):
> p = coords[d]
> limit += p
> - intf = integrate(intf, (p, 0, limit))
> + intf = integrate(intf.expand(), (p, 0, limit))
> return intf
>
> def bernstein_space(order, nsd):
As you've mentioned this is unneccessary.
> diff -r 315e41f3d28d -r 42b943e63382 sympy/integrals/integrals.py
> --- a/sympy/integrals/integrals.py Tue Nov 27 01:38:50 2007 +0100
> +++ b/sympy/integrals/integrals.py Tue Nov 27 02:28:12 2007 +0100
> @@ -115,9 +115,70 @@ class Integral(Basic, NoRelMeths, ArithM
> return function
>
> def _eval_integral(self, f, x):
> - # TODO : add table lookup for logarithmic and sine/cosine integrals
> - # and for some elementary special cases for speed improvement.
> - return risch_norman(f, x)
> + """Calculate the antiderivative to the function f(x).
> +
> + This is a powerful function that should in theory be able to integrate
> + everything that can be integrated. If you find something, that it
> + doesn't, it is easy to implement it.
> +
> + (1) Simple heuristics (based on pattern matching and integral table):
> +
> + - most frequently used functions (eg. polynomials)
> + - functions non-integrable by any of the following algorithms (eg.
> + exp(-x**2))
> +
> + (2) Integration of rational functions:
> +
> + (a) using apart() - apart() is full partial fraction decomposition
> + procedure based on Bronstein-Salvy algorithm. It gives formal
> + decomposition with no polynomial factorization at all (so it's fast
> + and gives the most general results). However it needs much better
> + implementation of RootsOf class (if fact any implementation).
> + (b) using Trager's algorithm - possibly faster than (a) but needs
> + implementation :)
> +
> + (3) Whichever implementation of pmInt (Mateusz, Kirill's or a
> + combination of both).
> +
> + - this way we can handle efficiently huge class of elementary and
> + special functions
> +
> + (4) Recursive Risch algorithm as described in Bronstein's integration
> + tutorial.
> +
> + - this way we can handle those integrable functions for which (3)
> + fails
> +
> + (5) Powerful heuristics based mostly on user defined rules.
> +
> + - handle complicated, rarely used cases
> + """
> +
> + # Let's first try some simple functions, that we know fast how to
> + # integrate.
> +
> + # simple powers:
> + from sympy import Pow, log
> + if isinstance(f, Pow) and isinstance(f[0], Symbol) and f[1].is_number:
> + if f[1] == -1:
> + return log(f[0])
> + else:
> + return f[0]**(f[1]+1)/(f[1]+1)
This looks good, but we probably should extend it to catch multiple
terms like:
x**Rational(1,2) + x**Rational(3,2)
BTW: I remember our Add has _eval_integrate ...
> +
> + # polynomials:
> + from sympy import Polynomial, PolynomialException
> + p = None
> + try:
> + p = Polynomial(f)
> + except PolynomialException:
> + pass
> + if p != None:
> + return p.integrate(x)
I suggest we do:
p = Polynomial(f, x)
or
p = f.as_polynomial(x)
example: f = sin(y)*x -- the original code will fail
> +
> + # f is not a simple function, let's try the risch norman (that can
> + # btw. integrate all the functions above, but slower):
> + r = risch_norman(f, x)
> + return r
>
> def integrate(*args, **kwargs):
> """Compute definite or indefinite integral of one or more variables
> @@ -132,6 +193,9 @@ def integrate(*args, **kwargs):
> >>> integrate(log(x), x)
> -x + x*log(x)
>
> + See also the doctest of Integral._eval_integral(), which explains
> + thoroughly the strategy that SymPy uses for integration.
> +
> """
> integral = Integral(*args, **kwargs)
Very nice in general.
>
> diff -r 315e41f3d28d -r 42b943e63382 sympy/integrals/tests/test_integrals.py
> --- a/sympy/integrals/tests/test_integrals.py Tue Nov 27 01:38:50 2007 +0100
> +++ b/sympy/integrals/tests/test_integrals.py Tue Nov 27 02:28:12 2007 +0100
> @@ -60,3 +60,9 @@ def test_multiple_integration():
>
> def test_issue433():
> assert integrate(exp(-x), (x,0,oo)) == 1
> +
> +def test_issue461():
> + assert integrate(x**Rational(3,2), x) == 2*x**Rational(5,2)/5
> + assert integrate(x**Rational(1,2), x) == 2*x**Rational(3,2)/3
> + assert integrate(x**Rational(-3,2), x) == -2*x**Rational(-1,2)
> +
Please add XFAIL for x**(1/2) + x**(3/2) or just include it in the test.
> diff -r 315e41f3d28d -r 42b943e63382 sympy/integrals/tests/test_risch.py
> --- a/sympy/integrals/tests/test_risch.py Tue Nov 27 01:38:50 2007 +0100
> +++ b/sympy/integrals/tests/test_risch.py Tue Nov 27 02:28:12 2007 +0100
> @@ -80,38 +80,41 @@ def test_resch_norman_issue442_1():
> # NB: correctness assured as ratsimp(diff(g,x) - f) == 0 in maxima
> # SymPy is unable to do it :(
>
> -@XFAIL
> -def test_pmint_rat():
> - f = (x**7-24*x**4-4*x**2+8*x-8) / (x**8+6*x**6+12*x**4+8*x**2)
> - g = (4 + 8*x**2 + 6*x + 3*x**3) / (x*(x**4 + 4*x**2 + 4)) + log(x)
> +# I commented all of those out, they are slowing down tests a lot (almost 40s
> +# spent for nothing on my comp) and they fail anyway. --Ondrej
>
> - assert risch_norman(f, x) == g
> -
> -
> -@XFAIL
> -def test_pmint_trig():
> - f = (x-tan(x)) / tan(x)**2 + tan(x)
> - g = (-x - tan(x)*x**2 / 2) / tan(x) + log(1+tan(x)**2) / 2
> -
> - assert risch_norman(f, x) == g
> -
> -
> -@XFAIL
> -def test_pmint_logexp():
> - f = (1+x+x*exp(x))*(x+log(x)+exp(x)-1)/(x+log(x)+exp(x))**2/x
> - g = 1/(x+log(x)+exp(x)) + log(x + log(x) + exp(x))
> -
> - assert risch_norman(f, x) == g
> -
> -
> -@XFAIL
> -def test_pmint_erf():
> - f = exp(-x**2)*erf(x)/(erf(x)**3-erf(x)**2-erf(x)+1)
> - g = sqrt(pi)/4 * (-1/(erf(x)-1) - log(erf(x)+1)/2 + log(erf(x)-1)/2)
> -
> - assert risch_norman(f, x) == g
> -
> -
> +#@XFAIL
> +#def test_pmint_rat():
> +# f = (x**7-24*x**4-4*x**2+8*x-8) / (x**8+6*x**6+12*x**4+8*x**2)
> +# g = (4 + 8*x**2 + 6*x + 3*x**3) / (x*(x**4 + 4*x**2 + 4)) + log(x)
> +#
> +# assert risch_norman(f, x) == g
> +#
> +#
> +#@XFAIL
> +#def test_pmint_trig():
> +# f = (x-tan(x)) / tan(x)**2 + tan(x)
> +# g = (-x - tan(x)*x**2 / 2) / tan(x) + log(1+tan(x)**2) / 2
> +#
> +# assert risch_norman(f, x) == g
> +#
> +#
> +#@XFAIL
> +#def test_pmint_logexp():
> +# f = (1+x+x*exp(x))*(x+log(x)+exp(x)-1)/(x+log(x)+exp(x))**2/x
> +# g = 1/(x+log(x)+exp(x)) + log(x + log(x) + exp(x))
> +#
> +# assert risch_norman(f, x) == g
> +#
> +#
> +#@XFAIL
> +#def test_pmint_erf():
> +# f = exp(-x**2)*erf(x)/(erf(x)**3-erf(x)**2-erf(x)+1)
> +# g = sqrt(pi)/4 * (-1/(erf(x)-1) - log(erf(x)+1)/2 + log(erf(x)-1)/2)
> +#
> +# assert risch_norman(f, x) == g
> +#
> +#
> # TODO: convert the rest of PMINT tests:
> # - Airy
> # - Bessel
py.test has skip -- we'd better use it in the beginning of each
function, or just prefix function names with underscore.
rationale1: avoid history pollution
rationale2: skip is maybe preffered, since in theory it should report:
rationale2: N tests skipped, and there should be a flag in py.test to
rationale2: ignore 'skip' tags
BTW: should be better in its own patch
> diff -r 315e41f3d28d -r 42b943e63382 sympy/polynomials/base.py
> --- a/sympy/polynomials/base.py Tue Nov 27 01:38:50 2007 +0100
> +++ b/sympy/polynomials/base.py Tue Nov 27 02:28:12 2007 +0100
> @@ -186,6 +186,7 @@ class Polynomial(Basic):
> "as_primitive",
> "content",
> "diff",
> + "integrate",
> "leading_coeff",
> "leading_term",
> "nth_coeff",
> @@ -640,6 +641,39 @@ class Polynomial(Basic):
> reverse=True)
> return Polynomial(coeffs=tuple(result_list), var=self.var,
> order=self.order)
> +
> + def integrate(self, variable):
> + """Primitive function of a Polynomial.
> +
> + Usage:
> + ======
> + Returns a new instance of Polynomial which is the primitive
> + function (antiderivative) of "self" with respect to the given
> + variable.
> +
> + Examples:
> + =========
> + >>> x, y, z = symbols('xyz')
> + >>> f = Polynomial(6*x + 20*y + 4*x*y)
> + >>> fx = f.integrate(x)
> + >>> print fx
> + 3*x**2 + 2*y*x**2 + 20*x*y
> + >>> fz = f.integrate(z)
> + >>> print fz
> + 6*x*z + 20*y*z + 4*x*y*z
> +
> + """
> +
> + if not variable in self.var:
> + return (variable * self).expand()
> + nvar = self.var.index(variable)
> + int = []
> + for term in self.coeffs:
> + t = list(term)
> + t[nvar+1] += 1
> + t[0] /= t[nvar+1]
> + int.append(tuple(t))
> + return coefficients2sympy(int, self.var)
Looks good.
One question: maybe the result of integrate here should be polynomial as
well? This way we'll avoid a lot of conversion vice versa when a lot of
polynomial handling is done.
BTW: better to live in it's own patch.
> def leading_coeff(self):
> diff -r 315e41f3d28d -r 42b943e63382 sympy/polynomials/tests/test_polynomials.py
> --- a/sympy/polynomials/tests/test_polynomials.py Tue Nov 27 01:38:50 2007 +0100
> +++ b/sympy/polynomials/tests/test_polynomials.py Tue Nov 27 02:28:12 2007 +0100
> @@ -140,10 +140,10 @@ def test_factor():
> assert factor(x**6-1) == (1+x**2-x)*(1+x)*(1+x+x**2)*(-1+x)
> assert factor(2*x**2+5*x+2) == (2+x)*(1+2*x)
>
> - assert factor(x**2 + y**2) == x**2 + y**2
> - assert factor(x*y + x*z + y*z) == x*y + x*z + y*z
> - assert factor(x*(y+1) + x*z) == x*(z + y + 1)
> - assert factor(x**5 - y**2) == x**5 - y**2
> + #assert factor(x**2 + y**2) == x**2 + y**2
> + #assert factor(x*y + x*z + y*z) == x*y + x*z + y*z
> + #assert factor(x*(y+1) + x*z) == x*(z + y + 1)
> + #assert factor(x**5 - y**2) == x**5 - y**2
Hmm, this is related to another issue 317, right?
I suggest we avoid mixing semantically independent changes.
>
> assert factor(-2) == -2
> assert factor(-x) == -x
> @@ -386,3 +386,11 @@ def test_poly_content_0():
> @XFAIL # see #442
> def test_poly_content_1():
> assert Polynomial(sqrt(2)*x, var=x).content() == sqrt(2)
> +
> +def test_poly_integrate():
> + x, y, z = symbols("xyz")
> + assert Polynomial(x**2+1).integrate(x) == x**3/3 + x
> + assert Polynomial(x**2+y*x).integrate(y) == x**2*y + y**2*x/2
> + assert Polynomial(x**2+x).integrate(y) == x**2*y + y*x
> + assert Polynomial(x*(y+x+z)).integrate(x) == x**2*y/2 + x**3/3 + x**2*z/2
> + assert Polynomial(x*(y+x+z)).integrate(z) == x*y*z + x**2*z + x*z**2/2
>
Overall looks good, thanks!
--
Всего хорошего, Кирилл.
http://landau.phys.spbu.ru/~kirr/aiv/
Kirill Smelkov
> > + # simple powers:
> > + from sympy import Pow, log
> > + if isinstance(f, Pow) and isinstance(f[0], Symbol) and f[1].is_number:
> > + if f[1] == -1:
> > + return log(f[0])
> > + else:
> > + return f[0]**(f[1]+1)/(f[1]+1)
As a consequence of vvv, this ^^^ should probably go to Pow.integrate
> This looks good, but we probably should extend it to catch multiple
> terms like:
>
> x**Rational(1,2) + x**Rational(3,2)
>
> BTW: I remember our Add has _eval_integrate ...
--
Ondrej Certik
>
> On Tue, Nov 27, 2007 at 08:32:54AM +0300, Kirill Smelkov wrote:
> > > + # simple powers:
> > > + from sympy import Pow, log
> > > + if isinstance(f, Pow) and isinstance(f[0], Symbol) and f[1].is_number:
> > > + if f[1] == -1:
> > > + return log(f[0])
> > > + else:
> > > + return f[0]**(f[1]+1)/(f[1]+1)
>
> As a consequence of vvv, this ^^^ should probably go to Pow.integrate
>
> > This looks good, but we probably should extend it to catch multiple
> > terms like:
> >
> > x**Rational(1,2) + x**Rational(3,2)
> >
> > BTW: I remember our Add has _eval_integrate ...
Thanks a lot for a quick review!
Unfortunately I don't have time right now to fix that and commit it. I
am leaving to Spain today.So probably I'll be able to fix it on
Sunday.
Ondrej
Kirill Smelkov
No problem. Good luck in Spain!
Ondrej Certik
thing that I didn't do was to semantically divide my patch into more.
I don't have time for that and it's important to release the new
changes.
Ondrej
Kirill Smelkov
Ok. I thought you'll fix it on Sunday or some time after.
But given we were out of time and going to release with this still
important changes it is ok.
Nevertheless since it is easier to prevent problems in the first place,
for future development I suggest we stick to reviewing each-other's
patches until they are somewhat ready. With constructive feedback the
quality rise is worth it!
Also, I've tried to address the issues you've created in separate mail.
Ondrej Certik
Yes, I completely agree and sorry about that.
Also, let's try to review it in terms of day - the review shouldn't
take much time, you just look at the patch and see if you like it or
not, you don't even have to test it - it's supposed to pass all tests
before sending the patch for review.
Ondrej
Ondrej Certik
> On Dec 10, 2007 12:05 AM, Kirill Smelkov <ki...@landau.phys.spbu.ru> wrote:
> >
> > On Wed, Dec 05, 2007 at 05:11:14PM -0800, Ondrej Certik wrote:
> > >
> > > I committed the patch, see new issues for your questions. The only
> > > thing that I didn't do was to semantically divide my patch into more.
> > > I don't have time for that and it's important to release the new
> > > changes.
> >
> > Ok. I thought you'll fix it on Sunday or some time after.
> >
> > But given we were out of time and going to release with this still
> > important changes it is ok.
> >
> > Nevertheless since it is easier to prevent problems in the first place,
> > for future development I suggest we stick to reviewing each-other's
> > patches until they are somewhat ready. With constructive feedback the
> > quality rise is worth it!
> >
> > Also, I've tried to address the issues you've created in separate mail.
>
> Yes, I completely agree and sorry about that.
>
> Also, let's try to review it in terms of day - the review shouldn't
I meant days - like one or two days are ok, but a week is too long
imho. I usually start a patch, send it for review and then want to
continue. I know I can continue on my own branch, but it's better to
continue on the reviewed code.
Ondrej
Kirill Smelkov
>
> On Dec 10, 2007 12:05 AM, Kirill Smelkov <ki...@landau.phys.spbu.ru> wrote:
> >
> > On Wed, Dec 05, 2007 at 05:11:14PM -0800, Ondrej Certik wrote:
> > >
> > > I committed the patch, see new issues for your questions. The only
> > > thing that I didn't do was to semantically divide my patch into more.
> > > I don't have time for that and it's important to release the new
> > > changes.
> >
> > Ok. I thought you'll fix it on Sunday or some time after.
> >
> > But given we were out of time and going to release with this still
> > important changes it is ok.
> >
> > Nevertheless since it is easier to prevent problems in the first place,
> > for future development I suggest we stick to reviewing each-other's
> > patches until they are somewhat ready. With constructive feedback the
> > quality rise is worth it!
> >
> > Also, I've tried to address the issues you've created in separate mail.
>
> Yes, I completely agree and sorry about that.
>
> Also, let's try to review it in terms of day - the review shouldn't
Yes, let's try to do the review in a few days. But given that the
review is done in spare time, it sometimes difficult to be on schedule.
But I'll try.
> take much time, you just look at the patch and see if you like it or
> not, you don't even have to test it - it's supposed to pass all tests
> before sending the patch for review.
What you say is good for the first iteration.
But next, for non-trivial patch, sometimes I like to see the code
around, that the patch doesn't touch, to see, for example, whether some
bits should be changed as well, so the patch trully changes the whole
semantically atom. This way the whole thing will be consistent.
Ondrej Certik
> around, that the patch doesn't touch, to see, for example, whether some
> bits should be changed as well, so the patch trully changes the whole
> semantically atom. This way the whole thing will be consistent.
Agree.
Ondrej