Skip to content

Commit f221aca

Browse files
mmaterammatera
committed
basic implementation of NIntegrate
1 parent 70bd781 commit f221aca

2 files changed

Lines changed: 298 additions & 1 deletion

File tree

‎mathics/builtin/exptrig.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -552,6 +552,7 @@ class ArcTanh(_MPMathFunction):
552552

553553
sympy_name = "atanh"
554554
mpmath_name = "atanh"
555+
numpy_name = "arctanh"
555556

556557
rules = {
557558
"Derivative[1][ArcTanh]": "1/(1-#^2)&",

‎mathics/builtin/numeric.py‎

Lines changed: 297 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -14,14 +14,15 @@
1414

1515
import sympy
1616
import mpmath
17+
import numpy as np
1718
from mpmath import mpf
1819
import math
1920
import hashlib
2021
import zlib
2122
import math
2223
from collections import namedtuple
2324
from contextlib import contextmanager
24-
from itertools import chain
25+
from itertools import chain, product
2526

2627

2728
from mathics.builtin.base import Builtin, Predefined
@@ -45,6 +46,7 @@
4546
Symbol,
4647
SymbolFalse,
4748
SymbolTrue,
49+
SymbolFailed,
4850
from_python,
4951
)
5052
from mathics.core.convert import from_sympy
@@ -255,6 +257,300 @@ def apply_other(self, expr, prec, evaluation):
255257
return Expression(head, *leaves)
256258

257259

260+
261+
def _internal_adaptative_simpsons_rule(f, a, b, **opts):
262+
"""
263+
1D adaptative Simpson's rule integrator
264+
Adapted from https://en.wikipedia.org/wiki/Adaptive_Simpson%27s_method
265+
by @mmatera
266+
267+
TODO: handle weak divergences
268+
"""
269+
wsr = 1./6.
270+
271+
tol = opts.get("tol")
272+
if not tol:
273+
tol = 1.e-10
274+
275+
maxrec = opts.get("maxrec")
276+
if not maxrec:
277+
maxrec = 150
278+
279+
def _quad_simpsons_mem(f, a, fa, b, fb):
280+
"""Evaluates the Simpson's Rule, also returning m and f(m) to reuse"""
281+
m = .5 * (a + b)
282+
try:
283+
fm = f(m)
284+
except ZeroDivisionError:
285+
fm = None
286+
287+
if fm is None or np.isinf(fm):
288+
m = m + 1e-10
289+
fm = f(m)
290+
return (m, fm, wsr * abs(b - a) * (fa + 4. * fm + fb))
291+
292+
def _quad_asr(f, a, fa, b, fb, eps, whole, m, fm, maxrec):
293+
"""
294+
Efficient recursive implementation of adaptive Simpson's rule.
295+
Function values at the start, middle, end of the intervals are retained.
296+
"""
297+
maxrec = maxrec - 1
298+
try:
299+
left = _quad_simpsons_mem(f, a, fa, m, fm)
300+
lm, flm, left = left
301+
right = _quad_simpsons_mem(f, m, fm, b, fb)
302+
rm, frm, right = right
303+
304+
delta = left + right - whole
305+
err = abs(delta)
306+
if err <= 15 * eps or maxrec == 0:
307+
return (left + right + delta / 15, err)
308+
left = _quad_asr(f, a, fa, m, fm, eps/2, left , lm, flm, maxrec)
309+
right = _quad_asr(f, m, fm, b, fb, eps/2, right, rm, frm, maxrec)
310+
return (left[0] + right[0], left[1] + right[1])
311+
except:
312+
raise
313+
314+
def ensure_evaluation(f, x):
315+
try:
316+
val = f(x)
317+
except ZeroDivisionError:
318+
return None
319+
if np.isinf(val):
320+
return None
321+
return val
322+
323+
invert_interval = False
324+
if a > b:
325+
b, a, invert_interval = a, b, True
326+
327+
fa, fb = ensure_evaluation(f, a), ensure_evaluation(f,b)
328+
if fa is None:
329+
x = 10. * machine_epsilon if a == 0 else a*(1. + 10.*machine_epsilon)
330+
fa = ensure_evaluation(f, x)
331+
if fa is None:
332+
raise Exception(f"Function undefined around {a}. Cannot integrate")
333+
if fb is None:
334+
x = -10. * machine_epsilon if b == 0 else b*(1. - 10.*machine_epsilon)
335+
fb = ensure_evaluation(f, x)
336+
if fb is None:
337+
raise Exception(f"Function undefined around {b}. Cannot integrate")
338+
339+
340+
m, fm, whole = _quad_simpsons_mem(f, a, fa, b, fb)
341+
return _quad_asr(f, a, fa, b, fb, tol, whole, m, fm, maxrec)
342+
343+
def _fubini(func, ranges, **opts):
344+
if not ranges:
345+
return 0.
346+
a, b = ranges[0]
347+
integrator = opts["integrator"]
348+
tol = opts.get("tol")
349+
if tol is None:
350+
opts["tol"] = 1.e-10
351+
tol = 1.e-10
352+
353+
if len(ranges) > 1:
354+
def subintegral(*u):
355+
def ff(*z):
356+
return func(*(u+z))
357+
358+
return _fubini(ff, ranges[1:], **opts)[0]
359+
opts["tol"] = 4. * tol
360+
return integrator(subintegral, a, b, **opts)
361+
else:
362+
return integrator(func, a, b, **opts)
363+
364+
365+
class NIntegrate(Builtin):
366+
"""
367+
<dl>
368+
<dt>'NIntegrate[$expr$, $interval$]'
369+
<dd>evaluates the definite integral of $expr$ numerically with a precision of $prec$ digits.
370+
</dl>
371+
372+
>> Table[1./NIntegrate[x^k,{x,0,1},Tolerance->1*^-6], {k,0,6}]
373+
= {1., 2., 3., 4., 5., 6., 7.}
374+
>> NIntegrate[Exp[-x],{x,0,Infinity},Tolerance->1*^-6]
375+
= 1.
376+
>> NIntegrate[Exp[x],{x,-Infinity, 0},Tolerance->1*^-6]
377+
= 1.
378+
>> NIntegrate[Exp[-x^2/2.],{x,-Infinity, Infinity},Tolerance->1*^-6]
379+
= 2.50663
380+
381+
Integrate singularities with weak divergences
382+
>> Table[NIntegrate[x^(1./k-1.),{x,0,1.},Tolerance->1*^-6], {k,1,7.}]
383+
= {1., 2., 3., 4., 5., 6., 7.}
384+
385+
Iterated Integral
386+
>> NIntegrate[x * y,{x, 0, 1},{y, 0, 1}]
387+
= 0.25
388+
"""
389+
# rules = {"N[Integrate[f_,{x_, a_Real, b__}]": "NIntegrate[f,{x,a,b}]",
390+
# "N[Integrate[f_,{x_, a_Real, b_Real}, prec_?NumericQ]":
391+
# 'NIntegrate[f,{x,a,b}, "AccuracyGoal"->prec]'
392+
# }
393+
394+
# messages = {'nnintegr':"The integrand `1` is not a numeric function",
395+
# 'nninflim': "The inferior limit `1` is not numeric",
396+
# 'nninflim': "The superior limit `1` is not numeric",
397+
# }
398+
399+
options = {'"Method"': '"Automatic"',
400+
'Tolerance': '1*^-10',
401+
'Accuracy': '1*^-10',
402+
'MaxRecursion': '10',
403+
}
404+
405+
406+
methods = {'Automatic': None,}
407+
408+
def __init__(self, *args, **kwargs):
409+
super().__init__(*args, **kwargs)
410+
self.methods["Internal"] = _internal_adaptative_simpsons_rule
411+
try:
412+
from scipy.integrate import (romberg, quad)
413+
print("Using scipy.integrate for numeric integration")
414+
self.methods["Quadrature"] = quad
415+
self.methods["Quadrature"] = romberg
416+
self.methods["Automatic"] = self.methods["Internal"]
417+
except Exception as e:
418+
print(e)
419+
print("Scipy not available. Using internal integrator")
420+
self.methods["Automatic"] = self.methods["Internal"]
421+
self.methods["Simpson"] = self.methods["Internal"]
422+
423+
@staticmethod
424+
def decompose_domain(interval, evaluation):
425+
if interval.has_form("System`Sequence", 1, None):
426+
intervals = []
427+
for leaf in interval.leaves:
428+
inner_interval = NIntegrate.decompose_domain(leaf, evaluation)
429+
if inner_interval:
430+
intervals.append(inner_interval)
431+
else:
432+
evaluation.message("Malformed interval")
433+
return None
434+
return intervals
435+
436+
if interval.has_form('System`List',3, None):
437+
intervals = []
438+
intvar = interval.leaves[0]
439+
if not isinstance(intvar, Symbol):
440+
return None
441+
boundaries = [a for a in interval.leaves[1:]]
442+
for i in range(len(boundaries)-1):
443+
intervals.append((boundaries[i],boundaries[i+1]))
444+
if len(intervals)>0:
445+
return (intvar, intervals)
446+
447+
return None
448+
evaluation.message("Domain not properly formatted")
449+
return None
450+
451+
def apply_2(self, func, domain, evaluation, options):
452+
"NIntegrate[func_, domain__, OptionsPattern[%(name)s]]"
453+
method = options['System`"Method"'].value
454+
tolerance = options['System`Tolerance'].evaluate(evaluation)
455+
tolerance = float(tolerance.value)
456+
accuracy = options['System`Accuracy'].evaluate(evaluation)
457+
accuracy = accuracy.value
458+
maxrecursion = options['System`MaxRecursion'].evaluate(evaluation)
459+
maxrecursion = maxrecursion.value
460+
461+
nintegrate_method = self.methods.get(method)
462+
if nintegrate_method is None:
463+
nintegrate_method = self.methods.get("Automatic")
464+
465+
domain = self.decompose_domain(domain, evaluation)
466+
if not domain:
467+
return
468+
if not isinstance(domain,list):
469+
domain = [domain]
470+
471+
coords = [axis[0] for axis in domain]
472+
intvars = Expression("List", *coords)
473+
integrand = Expression("Compile", intvars, func).evaluate(evaluation)
474+
if len(integrand.leaves)>=3:
475+
integrand = integrand.leaves[2].cfunc
476+
else:
477+
print(func, " could not be compiled")
478+
return
479+
results = []
480+
for subdomain in product(*[axis[1] for axis in domain]):
481+
# On each subdomain, check if the region is bounded.
482+
# If not, implement a coordinate map
483+
func2 = integrand
484+
subdomain2 = []
485+
coordtransform = []
486+
nulldomain = False
487+
for i, r in enumerate(subdomain):
488+
a = r[0].evaluate(evaluation)
489+
b = r[1].evaluate(evaluation)
490+
if a == b:
491+
nulldomain = True
492+
break
493+
elif a.get_head_name() == 'System`DirectedInfinity':
494+
if b.get_head_name() == 'System`DirectedInfinity':
495+
a = a.to_python()
496+
b = b.to_python()
497+
l = 1 - machine_epsilon
498+
if a == b:
499+
nulldomain = True
500+
break
501+
elif a < b:
502+
subdomain2.append([-l, l])
503+
else:
504+
subdomain2.append([l, -l])
505+
coordtransform.append((np.arctanh, lambda u:1./(1.-u**2)))
506+
else:
507+
if not b.is_numeric():
508+
evaluation.message("Boundary is not numeric.")
509+
return SymbolFailed
510+
z = a.leaves[0].value
511+
b = b.value
512+
subdomain2.append( [machine_epsilon, 1.])
513+
coordtransform.append((lambda u: b-z + z/u, lambda u: -z*u**(-2.)))
514+
elif b.get_head_name() == 'System`DirectedInfinity':
515+
if not a.is_numeric():
516+
evaluation.message("Boundary is not numeric.")
517+
return SymbolFailed
518+
a = a.value
519+
z = b.leaves[0].value
520+
subdomain2.append([machine_epsilon, 1.])
521+
coordtransform.append((lambda u: a-z + z/u, lambda u: z*u**(-2.)))
522+
elif a.is_numeric() and b.is_numeric():
523+
a = Expression("N", a).evaluate(evaluation).value
524+
b = Expression("N", b).evaluate(evaluation).value
525+
subdomain2.append([a, b])
526+
coordtransform.append(None)
527+
else:
528+
evaluation.message("Boundary is not numeric.")
529+
return SymbolFailed
530+
531+
if nulldomain:
532+
continue
533+
534+
if any(coordtransform):
535+
func2 = lambda *u: (integrand(*[x[0](u[i]) if x else u[i]
536+
for i, x in enumerate(coordtransform)])
537+
* np.prod([ jac[1](u[i] )
538+
for i, jac in enumerate(coordtransform) if jac]))
539+
opts = {#"" : accuracy,
540+
"tol" : tolerance,
541+
#"maxrec" : maxrecursion,
542+
}
543+
if len(subdomain2) >1:
544+
val = _fubini(func2, subdomain2,
545+
integrator=nintegrate_method, **opts)
546+
else:
547+
val = nintegrate_method(func2, *(subdomain2[0]), **opts)
548+
results.append(val)
549+
result = sum([r[0] for r in results])
550+
error = sum([r[1] for r in results])
551+
# print("result=", result, " error=",error)
552+
return from_python(result)
553+
258554
class MachinePrecision(Predefined):
259555
"""
260556
<dl>

0 commit comments

Comments
 (0)