Skip to content

Commit a37ffb2

Browse files
author
Release Manager
committed
gh-39094: Fix numerical_integral with parameters Fixes #39088. As pointed out by @neldredge-unco in #39088, supplying a value for the keyword argument `params` of `numerical_integral` yields a bunch of blank lines, and then a value of `0.0` for the integral. **Diagnosis.** The value of the integrand at each point is calculated by the function `c_f` (lines 49-62 of `src/sage/calculus/integration.pyx`): ``` cdef double c_f(double t, void *params) noexcept: cdef double value cdef PyFunctionWrapper wrapper wrapper = <PyFunctionWrapper> params try: if len(wrapper.the_parameters) != 0: value = wrapper.the_function(t, wrapper.the_parameters) else: value = wrapper.the_function(t) except Exception as msg: print(msg) return 0 return value ``` If an exception is raised while trying to calculate the value, then the error message is printed, and 0 is returned as the value of the function. Since parameters are not handled properly by the `numerical_integral` function, an exception is raised every time, but the error message is empty, so we get lots of blank lines and the integral is 0. **Remedy implemented in this PR.** 1. If `msg` has only whitespace, the PR instead prints an actual error message. This replaces the blank lines with something more helpful. 2. The function `c_f` is tagged `noexcept`, so we cannot allow an exception to be raised while printing the error message. Therefore, the PR puts the printing in a `try...except` block to catch all exceptions. 3. `numerical_integral` converts the integrand to a `fast_callable` function with only one variable (line 322 of `src/sage/calculus/integration.pyx`): ``` func = fast_callable(func, vars=[v], domain=float) ``` However, `numerical_integral` still passes `params` (the parameters) to the function `c_f`, so `c_f` is receiving more than one input value for a one-variable function. This causes an exception every time. To fix this, the PR executes `params = []` immediately after `func` is made into a `fast_callable`. This means that no unwanted extra parameters will be sent to `c_f`. 4. The doctest for `numerical_integral` was marked `# random` so it was essentially ignored during doctesting. (I think that is why this bug was not discovered long ago.) The PR changes this doctest to use a tolerance, so incorrect behaviour will be detected in the future. **Note.** The function `c_f` uses `except Exception`, instead of only catching a specific list of exceptions. This is not good python, but I think fixing this would belong on a different ticket, because it will probably not be trivial (since a `noexcept` function cannot raise exceptions). ### 📝 Checklist <!-- Put an `x` in all the boxes that apply. --> - [x] The title is concise and informative. - [x] The description explains in detail what this PR is about. - [x] I have linked a relevant issue or discussion. - [x] I have created tests covering the changes. - [x] I have updated the documentation and checked the documentation preview. ### ⌛ Dependencies <!-- List all open PRs that this PR logically depends on. For example, --> <!-- - #12345: short description why this is a dependency --> <!-- - #34567: ... --> URL: #39094 Reported by: DaveWitteMorris Reviewer(s): Travis Scrimshaw
2 parents 2059f42 + fa767fd commit a37ffb2

File tree

1 file changed

+23
-12
lines changed

1 file changed

+23
-12
lines changed

src/sage/calculus/integration.pyx

Lines changed: 23 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -59,7 +59,13 @@ cdef double c_f(double t, void *params) noexcept:
5959
else:
6060
value = wrapper.the_function(t)
6161
except Exception as msg:
62-
print(msg)
62+
try:
63+
if str(msg).strip():
64+
print(msg)
65+
else:
66+
print(f"Unable to evaluate function at {t}")
67+
except Exception:
68+
pass
6369
return 0
6470

6571
return value
@@ -144,17 +150,19 @@ def numerical_integral(func, a, b=None,
144150
For a Python function with parameters::
145151
146152
sage: f(x,a) = 1/(a+x^2)
147-
sage: [numerical_integral(f, 1, 2, max_points=100, params=[n]) for n in range(10)] # random output (architecture and os dependent)
148-
[(0.49999999999998657, 5.5511151231256336e-15),
149-
(0.32175055439664557, 3.5721487367706477e-15),
150-
(0.24030098317249229, 2.6678768435816325e-15),
151-
(0.19253082576711697, 2.1375215571674764e-15),
152-
(0.16087527719832367, 1.7860743683853337e-15),
153-
(0.13827545676349412, 1.5351659583939151e-15),
154-
(0.12129975935702741, 1.3466978571966261e-15),
155-
(0.10806674191683065, 1.1997818507228991e-15),
156-
(0.09745444625548845, 1.0819617008493815e-15),
157-
(0.088750683050217577, 9.8533051773561173e-16)]
153+
sage: [numerical_integral(f, 1, 2, max_points=100, params=[n])[0] # abs tol 1.0e-6
154+
....: for n in range(10)]
155+
[0.5000000000000000,
156+
0.3217505543966422,
157+
0.24030098317248832,
158+
0.19253082576711372,
159+
0.1608752771983211,
160+
0.138275456763492,
161+
0.1212997593570257,
162+
0.10806674191683492,
163+
0.09745444625553161,
164+
0.08875068305030848]
165+
158166
sage: y = var('y')
159167
sage: numerical_integral(x*y, 0, 1)
160168
Traceback (most recent call last):
@@ -323,6 +331,9 @@ def numerical_integral(func, a, b=None,
323331
if ell.is_numeric() and not ell.is_zero():
324332
raise ValueError('integral does not converge at infinity')
325333
func = fast_callable(func, vars=[v], domain=float)
334+
# `func` is now a function of one variable,
335+
# so it no longer needs any parameters
336+
params = []
326337

327338
if not isinstance(func, compiled_integrand):
328339
wrapper = PyFunctionWrapper()

0 commit comments

Comments
 (0)