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
-~----------~----~----~----~------~----~------~--~---
