Ondrej Certik wrote: > if the method='scipy' is specified, the scipy.integrate.quad function is > called > on the lambdified argument of the integral. > > Tests were added to sympy/test_external/test_scipy.py > > Signed-off-by: Ondrej Certik <ond...@certik.cz> > --- > sympy/core/evalf.py | 43 +++++++++++++++++++++------- > sympy/core/numbers.py | 2 +- > sympy/integrals/tests/test_integrals.py | 6 ++++ > sympy/test_external/test_scipy.py | 46 > +++++++++++++++++++++++++++++++ > 4 files changed, 85 insertions(+), 12 deletions(-) > create mode 100644 sympy/test_external/test_scipy.py > > diff --git a/sympy/core/evalf.py b/sympy/core/evalf.py > index a32d9b8..0314a03 100644 > --- a/sympy/core/evalf.py > +++ b/sympy/core/evalf.py > @@ -280,7 +280,6 @@ def add_terms(terms, prec, target_prec): > sum_accuracy = sum_exp + sum_bc - absolute_error > r = normalize(sum_sign, sum_man, sum_exp, sum_bc, target_prec, > round_nearest), sum_accuracy > - #print "returning", to_str(r[0],50), r[1] > return r > > def evalf_add(v, prec, options): > @@ -730,16 +729,34 @@ def f(t): > return result > > def evalf_integral(expr, prec, options): > - workprec = prec > - i = 0 > - maxprec = options.get('maxprec', INF) > - while 1: > - result = do_integral(expr, workprec, options) > - accuracy = complex_accuracy(result) > - if accuracy >= prec or workprec >= maxprec: > - return result > - workprec += prec - max(-2**i, accuracy) > - i += 1 > + method = options.get('method', 'mpmath') > + if method == 'scipy': > + from scipy.integrate import quad > + func = expr.args[0] > + x, (xlow, xhigh) = expr.args[1][0]
this will fail silently (and give a wrong result) for integration on several variables: >>> i = Integral(1, (x, 0, 1), (y, 0, oo)) >>> i In [32]: i Out[32]: oo 1 / / | | | | 1 dx dy | | / / 0 0 >>> i.evalf() >>> 1.00000000000000 Also, I think that x, (xlow, xhigh) = expr.limits[0] would be more elegant > + from sympy import lambdify > + f = lambdify(x, func, ["math"]) > + # scipy only operates on the standard floating point precision: > + eps = 1e-9 > + re, im = quad(f, xlow, xhigh, epsabs=eps, epsrel=eps) > + # chop imaginary parts that are below the precision > + if abs(im) < eps: > + im = 0 > + re = from_float(re) > + im = from_float(im) > + return re, im, re[3], im[3] > + elif method == 'mpmath': > + workprec = prec > + i = 0 > + maxprec = options.get('maxprec', INF) > + while 1: > + result = do_integral(expr, workprec, options) > + accuracy = complex_accuracy(result) > + if accuracy >= prec or workprec >= maxprec: > + return result > + workprec += prec - max(-2**i, accuracy) > + i += 1 > + raise ValueError("Unknown method: %s" % method) > > def check_convergence(numer, denom, n): > """ > @@ -1003,6 +1020,10 @@ def Basic_evalf(x, n=15, **options): > verbose=<bool> > Print debug information (default=False) > > + method=<str> > + mpmath ... uses mpmath for evaluation > + scipy .... uses scipy where possible (e.g. evaluation of > integrals) > + > """ > if not evalf_table: > _create_evalf_table() > diff --git a/sympy/core/numbers.py b/sympy/core/numbers.py > index f9cc5f2..b6a74c6 100644 > --- a/sympy/core/numbers.py > +++ b/sympy/core/numbers.py > @@ -642,7 +642,7 @@ def __ne__(self, other): > if other.is_comparable and not isinstance(other, Rational): other = > other.evalf() > if isinstance(other, Number): > if isinstance(other, Real): > - return bool(not mlib.feq(self._as_mpf_val(other._prec), > other._mpf_)) > + return bool(not mlib.mpf_eq(self._as_mpf_val(other._prec), > other._mpf_)) > return bool(self.p!=other.p or self.q!=other.q) > > return True # Rational != non-Number > diff --git a/sympy/integrals/tests/test_integrals.py > b/sympy/integrals/tests/test_integrals.py > index 301090c..9a7bc63 100644 > --- a/sympy/integrals/tests/test_integrals.py > +++ b/sympy/integrals/tests/test_integrals.py > @@ -339,3 +339,9 @@ def test_subs5(): > def test_integration_variable(): > raises(ValueError, "Integral(exp(-x**2), 3)") > raises(ValueError, "Integral(exp(-x**2), (3, -oo, oo))") > + > +def test_integrals_evalf(): > + i = Integral(exp(-x**2), (x, 0, oo)) > + assert NS(i.evalf()) == "0.886226925452758" > + assert NS(i.evalf(method="mpmath")) == "0.886226925452758" > + raises(ValueError, "i.evalf(method='nonsense5')") > diff --git a/sympy/test_external/test_scipy.py > b/sympy/test_external/test_scipy.py > new file mode 100644 > index 0000000..3525f9b > --- /dev/null > +++ b/sympy/test_external/test_scipy.py > @@ -0,0 +1,46 @@ > +# This testfile tests SymPy <-> SciPy compatibility > + > +# Don't test any SymPy features here. Just pure interaction with SymPy. > +# Always write regular SymPy tests for anything, that can be tested in pure > +# Python (without numpy). Here we test everything, that a user may need when > +# using SymPy with NumPy > + > +try: > + import scipy > +except ImportError: > + #py.test will not execute any tests now > + disabled = True > + > + > +from sympy import (Integral, exp, oo, Symbol, sstr, sympify, raises) > + > +def NS(e, n=15, **options): > + return sstr(sympify(e).evalf(n, **options), full_prec=True) not really obvious what this function does, a docstring would be appreciated > + > +def test_integrals_evalf_scipy1(): > + x = Symbol("x") > + i = Integral(exp(-x**2), (x, 0, oo)) > + assert NS(i.evalf()) == "0.886226925452758" > + assert NS(i.evalf(method="mpmath")) == "0.886226925452758" > + assert NS(i.evalf(method="scipy")) == "0.886226925452758" > + > +def test_integrals_evalf_scipy2(): > + x = Symbol("x") > + i = Integral(exp(-x**2), (x, -oo, 0)) > + assert NS(i.evalf()) == "0.886226925452758" > + assert NS(i.evalf(method="mpmath")) == "0.886226925452758" > + assert NS(i.evalf(method="scipy")) == "0.886226925452758" > + > +def test_integrals_evalf_scipy3(): > + x = Symbol("x") > + i = Integral(exp(-x**2), (x, -oo, oo)) > + assert NS(i.evalf()) == "1.77245385090552" > + assert NS(i.evalf(method="mpmath")) == "1.77245385090552" > + assert NS(i.evalf(method="scipy")) == "1.77245385090552" > + > +def test_integrals_evalf_scipy4(): > + x = Symbol("x") > + i = Integral(exp(-x**2), (x, -0.1, 0.1)) > + assert NS(i.evalf()) == "0.199335328580673" > + assert NS(i.evalf(method="mpmath")) == "0.199335328580673" > + assert NS(i.evalf(method="scipy")) == "0.199335328580673" --~--~---------~--~----~------------~-------~--~----~ You received this message because you are subscribed to the Google Groups "sympy-patches" group. To post to this group, send email to sympy-patches@googlegroups.com To unsubscribe from this group, send email to sympy-patches+unsubscr...@googlegroups.com For more options, visit this group at http://groups.google.com/group/sympy-patches?hl=en -~----------~----~----~----~------~----~------~--~---