Object-Oriented Programming¶
In this notebook, I attempt to explain by example the utility of object oriented programming. In particular, I highlight the use of inheritance to ensure a common interface to objects. This concept is illustrating by focusing on numerical integration. We have by now a large number of different algorithms to numerically approximate integrals. It makes sense that we should be able to access these algorithms in similar ways, and using that similarity and thinking abstractly is the key to being able to build up structure to solve more advanced problems.
This discussion is somewhat similar to Chapter 9 of
- [L] A Primer on Scientific Programming with Python by Hans Petter Langtangen, 2nd edition.
This Chapter has further examples and goes into greater depth, while this notebook focuses on one family of examples.
# Standard imports
import numpy as np
import matplotlib.pyplot as plt
import math as m
from mpmath import mp, iv
from numpy import random
import random as random_number
Inheritance and integration methods¶
We have several integration approximation methods. Here is one, implemented as a class:
class MidpointRule:
def __init__(self, a, b):
self._a = a
self._b = b
def integrate(self, f):
a = self._a
b = self._b
return (b-a) * f((a+b)/2)
Here we demonstrate it's use to approximate $$\int_{0}^\pi \sin x~dx=2.$$
mr = MidpointRule(0.0, m.pi)
f = lambda x: m.sin(x)
mr.integrate(f)
Here is another rule we learned about, implemented in almost the same way:
class TrapezoidRule:
def __init__(self, a, b):
self._a = a
self._b = b
def integrate(self, f):
a = self._a
b = self._b
return (b-a)/2 * (f(a)+f(b))
tr = TrapezoidRule(0.0, m.pi)
f = lambda x: m.sin(x)
tr.integrate(f)
Both the midpoint rule and trapezoid rule are examples of integration approximation rules. We have other examples too: Simpson's rule, Gaussian quadrature, composite rules, and most recently Monte Carlo Integration. Because the are all integration approximation rules, they all integrate functions, so in their class implementations we would always have an integrate(f) method taking a function as input.
It is useful to program all these class in a similar way (like above) so that they can be used interchangably. This interchangability can be guaranteed through inheritance, which also offers other benefits including code reuse. For example, the handling of endpoints is the same in the above two classes.
Here is an abstract base class representing all the common features we might expect from an integration method for approximating integrals in an interval $[a,b]$.
class IntervalIntegrationMethod:
def __init__(self, a, b):
"""Represent an integral over [a, b]."""
self._a = a
self._b = b
def left_endpoint(self):
"""Return the left endpoint of the interval."""
return self._a
def right_endpoint(self):
"""Return the right endpoint of the interval."""
return self._b
def interval(self):
"""Return a pair consisting of the left and right endpoints."""
return (self._a, self._b)
# Abstract Method
def integrate(self, f):
"""Return the approximation of the integral of f."""
raise NotImplementedError("integrate(f) was not implemented")
Note that the __init__ method manages the endpoints. The three methods left_endpoint(), right_endpoint() and interval() look normal and will work fine. The integrate(f) method will raise an error if called.
The class doesn't do much, but we can do this for example:
iim = IntervalIntegrationMethod(1, 2)
iim.interval()
On the other hand, calling
iim.integrate(lambda x: x**2)
Will lead to a NotImplementedError.
We will now recreate a MidpointRule class that inherits from IntervalIntegrationMethod and implement the integrate(f) method.
class MidpointRule(IntervalIntegrationMethod):
def integrate(self, f):
"""Return the approximation of the integral of f."""
a = self._a
b = self._b
return (b-a) * f((a+b)/2)
The indication that we inherit from IntervalIntegrationMethod is in the class declaration:
class MidpointRule(IntervalIntegrationMethod):
By including the integrate(f) method in the definition, we override the previous definition (which led purposely to errors). When we call integrate(f) from a MidpointRule object, the new method will be called.
The class can be used just like before.
mr = MidpointRule(0.0, m.pi)
f = lambda x: m.sin(x)
mr.integrate(f)
The other common rules can be implemented just as quickly:
class TrapezoidRule(IntervalIntegrationMethod):
def integrate(self, f):
"""Return the approximation of the integral of f."""
a = self._a
b = self._b
return (b-a)/2 * (f(a)+f(b))
class SimpsonsRule(IntervalIntegrationMethod):
def integrate(self, f):
"""Return the approximation of the integral of f."""
a = self._a
b = self._b
return (b-a)/6 * (f(a)+4*f((a+b)/2)+f(b))
Simpson's rule actually gives a fairly close answer to our test function.
mr = SimpsonsRule(0.0, m.pi)
f = lambda x: m.sin(x)
mr.integrate(f)
Recall that a composite rule is formed by subdividing an interval $[a,b]$ into $n$ subintervals, typically of equal size. The endpoints can be found via
np.linspace(a, b, n+1)
For example:
np.linspace(0, 1, 10+1)
We can create a new composite integration rule for integrals over [a,b] with the following data: a, b, a number of integrals n and a basic integration rule as above. We can accept this data by overriding the __init__ method. But, the new __init__ method needs to call the original __init__ method to handle the endpoints. This can be done and here is an example:
class CompositeRule(IntervalIntegrationMethod):
def __init__(self, a, b, n, BasicIntegrationMethod):
"""Construct a new CompositeRule for integrals over [a, b] by dividing
the interval into n equal sized subintervals and applying the
BasicIntegrationMethod on each one.
BasicIntegrationMethod(x,y) should return an integration method over
the interval [x, y].
"""
# Call the __init__ method in IntervalIntegrationMethod
# to set up the endpoints:
IntervalIntegrationMethod.__init__(self, a, b)
points = np.linspace(a, b, n+1)
self._i = [] # Will store a list of integrators over subintervals
for j in range(n):
self._i.append( BasicIntegrationMethod(points[j], points[j+1]) )
def integrate(self, f):
"""Return the approximation of the integral of f."""
total = 0.0
for integrator in self._i:
total += integrator.integrate(f)
return total
cr = CompositeRule(0.0, m.pi, 10, MidpointRule)
f = lambda x: m.sin(x)
cr.integrate(f)
That faired only a little better than Simpson's rule. Because of the flexibility of our CompositeRule class, we can also create a composite Simpson's rule:
cr = CompositeRule(0.0, m.pi, 10, SimpsonsRule)
f = lambda x: m.sin(x)
cr.integrate(f)
As another example, we can implement the Monte Carlo integration method:
class MonteCarloInterval(IntervalIntegrationMethod):
def __init__(self, a, b, trials):
"""The MonteCarlo Integration method for integrating functions in
[a, b], using `trials` points taken at random."""
# Call the __init__ method in IntervalIntegrationMethod
# to set up the endpoints:
IntervalIntegrationMethod.__init__(self, a, b)
self._trials = trials
def integrate(self, f):
"""Return the approximation of the integral of f."""
a = self._a
b = self._b
points = (b-a)*np.random.random_sample(self._trials) + a
return (b-a) * np.sum(f(points)) / self._trials
We print out $10$ approximates for $\int_0^\pi \sin x~dx$ belo using 100 sample points.
mci = MonteCarloInterval(0.0, np.pi, 100)
f = lambda x: np.sin(x)
for i in range(10):
print(mci.integrate(f))
We can even apply the composite method to our Monte Carlo integration implementation. We will divide $[0, \pi]$ into $10$ intervals and then apply $10$ point Monte Carlo integration to each subinterval. This will still use $10$ points but they will be more evenly distributed, so we expect better results.
Note that our class note that our class MonteCarloInterval takes an extra parameter to initialize an object, the number of points. The CompositeRule expects to be passed a rule which only takes the endpoints as inputs. We can substute a function for this rule. See the example below.
def ten_point_mci(a, b):
"""Return the MonteCarloInterval integration method on [a, b] using 10 points."""
return MonteCarloInterval(a, b, 10)
cr = CompositeRule(0.0, np.pi, 10, ten_point_mci)
f = lambda x: np.sin(x)
for i in range(10):
print(cr.integrate(f))
Abstracting Further: Iterated integrals¶
It would be nice not to be restricted to intervals in the real line. We could for example hope to integrate functions over regions in the plane or in $3$-dimensional space.
We'd want an IntegrationMethod to again have an integrate(f) method, but it would be nice to know what $f$ takes as input. It could be for example be a real number or a vector. As far as this notebook is concerned, we will fix the domain to be ${\mathbb R}^n$ for some $n$, so we will add a new method domain_dimension() which returns this dimension $n$.
Below is our abstract IntegrationMethod class.
class IntegrationMethod:
# Abstract Method
def domain_dimension(self):
"""Return the dimension of the domain of functions being integrated."""
raise NotImplementedError("domain_dimension() was not implemented")
# Abstract Method
def integrate(self, f):
"""Return the approximation of the integral of f."""
raise NotImplementedError("integrate(f) was not implemented")
We have rewritten our IntervalIntegrationMethod to inherit from this class. The changes are listed below:
- The class is now declared
class IntervalIntegrationMethod(IntegrationMethod):
- We no longer need to include the
integrate(f)method since it is already inIntegrationMethod. - We implement the
domain_dimension()method, which just returns $1$.
class IntervalIntegrationMethod(IntegrationMethod):
def __init__(self, a, b):
"""Represent an integral over [a, b]."""
self._a = a
self._b = b
def left_endpoint(self):
"""Return the left endpoint of the interval."""
return self._a
def right_endpoint(self):
"""Return the right endpoint of the interval."""
return self._b
def interval(self):
"""Return a pair consisting of the left and right endpoints."""
return (self._a, self._b)
def domain_dimension(self):
"""Return the dimension of the domain of functions being integrated."""
return 1
Below I've included the various classes we have already defined. They have not changed, but we need to rerun the code so that they inherit from the new classes above.
class MidpointRule(IntervalIntegrationMethod):
def integrate(self, f):
"""Return the approximation of the integral of f."""
a = self._a
b = self._b
return (b-a) * f((a+b)/2)
class TrapezoidRule(IntervalIntegrationMethod):
def integrate(self, f):
"""Return the approximation of the integral of f."""
a = self._a
b = self._b
return (b-a)/2 * (f(a)+f(b))
class SimpsonsRule(IntervalIntegrationMethod):
def integrate(self, f):
"""Return the approximation of the integral of f."""
a = self._a
b = self._b
return (b-a)/6 * (f(a)+4*f((a+b)/2)+f(b))
class CompositeRule(IntervalIntegrationMethod):
def __init__(self, a, b, n, BasicIntegrationMethod):
"""Construct a new CompositeRule for integrals over [a, b] by dividing
the interval into n equal sized subintervals and applying the
BasicIntegrationMethod on each one.
BasicIntegrationMethod(x,y) should return an integration method over
the interval [x, y].
"""
# Call the __init__ method in IntervalIntegrationMethod
# to set up the endpoints:
IntervalIntegrationMethod.__init__(self, a, b)
points = np.linspace(a, b, n+1)
self._i = [] # Will store a list of integrators over subintervals
for j in range(n):
self._i.append( BasicIntegrationMethod(points[j], points[j+1]) )
def integrate(self, f):
"""Return the approximation of the integral of f."""
total = 0.0
for integrator in self._i:
total += integrator.integrate(f)
return total
class MonteCarloInterval(IntervalIntegrationMethod):
def __init__(self, a, b, trials):
"""The MonteCarlo Integration method for integrating functions in
[a, b], using `trials` points taken at random."""
# Call the __init__ method in IntervalIntegrationMethod
# to set up the endpoints:
IntervalIntegrationMethod.__init__(self, a, b)
self._trials = trials
def integrate(self, f):
"""Return the approximation of the integral of f."""
a = self._a
b = self._b
points = (b-a)*np.random.random_sample(self._trials) + a
return (b-a) * np.sum(f(points)) / self._trials
The main manner in which we extend integration to higher dimensions is by iterated integration. For example, to integrate over a rectangle $R=[a,b] \times [c,d]$ we can do either: $$\int_R f~dA= \int_a^b \int_c^d f(x,y)~dy~dx = \int_c^d \int_a^b f(x,y)~dx~dy.$$ (This is true for $f$ continuous by Fubini's theorem.)
Thinking more abstractly, we might have a region $R$ in ${\mathbb R}^m$ and a region $S$ in ${\mathbb R}^n$. Then the product $R \times S$ lives in ${\mathbb R}^{m+n}$ and consists of vectors of the form $(v_1, v_2, \ldots, v_m, w_1, w_2, \ldots, w_n)$ with $(v_1, v_2, \ldots, v_m) \in R$ and $(w_1, \ldots, w_n) \in S$. We have $$\int_{R \times S} f(v,w) dVol = \int_R \int_S f(v,w)~dw~dv = \int_S \int_R f(v,w)~dv~dw.$$ Here $dV_{m}$ means we integrate with respect to volume in $m$-dimensional space.
Remark. This is abstract, but you can think about it in examples. If $R$ is the unit disk in the plane and $S$ is an interval, then $R \times S$ is a solid cylinder. It is useful to think in generality, because it avoids duplication of code.
To prepare, we need to think about the issue of what the integral $$\int_R \int_S f(v,w)~dw~dv$$ means and how to code it. The outer integral fixes a $v$ value. Then the inner integral is really computing the integral of the function $$g_v(w) = f(v,w)$$ with respect to $w$. Then we can think of the inner integral as a function of $v$: $$h(v) = \int_S g_v(w)~dw$$ Then we want to integrate $h(v)$ with respect to $v$ and this is the meaning of the double integral.
v = np.array([1, 1])
w = np.array([2, 2, 2])
np.concatenate([v, w])
np.array([1])
np.concatenate((np.array([1]),np.array([1])))
np.append(v,w)
class IteratedIntegral(IntegrationMethod):
def __init__(self, outer_integral, inner_integral):
"""Construct an iterated integral from two integration methods.
When we integrate a function, the outer integral will integrate over
the first coordinates, and the inner integral will iterate over the remaining
coordinates.
"""
self._outer = outer_integral
self._inner = inner_integral
def domain_dimension(self):
"""Return the dimension of the domain of functions being integrated."""
return self._outer.domain_dimension() + self._inner.domain_dimension()
def integrate(self, f):
"""Return the approximation of the integral of f."""
def h(v):
def g(w):
return f(np.append(v,w))
return self._inner.integrate(g)
return self._outer.integrate(h)
Test 1. To test our class we will verify that $$\iint_R x y^2 ~dA = \frac{2}{3}$$ when $R$ is the rectangle defined by $0 \leq x \leq 2$ and $0 \leq y \leq 1$. As an iterated integral, this is $$\int_0^2 \int_0^1 x y^2 ~dy~dx = \frac{2}{3}.$$
For the inner and outer methods, we will use composite trapezoid rules with $10$ subintervals.
outer_integral_method = CompositeRule(0, 2, 10, TrapezoidRule)
inner_integral_method = CompositeRule(0, 1, 10, TrapezoidRule)
integral_method = IteratedIntegral(outer_integral_method,
inner_integral_method)
Our function will take a numpy array with two coordinates as its input.
def f(v):
x,y = v
return x * y**2
integral_method.integrate(f)
Remark. The above method amounts to a quadrature rule for the $[0,2] \times [0,1]$ rectangle using an $11 \times 11$ grid of points.
If we use Simpson's rule instead, we should get an exact answer (since the degree of accuracy of Simpson's rule is $3$).
outer_integral_method = CompositeRule(0, 2, 10, SimpsonsRule)
inner_integral_method = CompositeRule(0, 1, 10, SimpsonsRule)
integral_method = IteratedIntegral(outer_integral_method,
inner_integral_method)
def f(v):
x,y = v
return x * y**2
integral_method.integrate(f)
Test 2. Now we will demonstrate how we can compute a triple integral. Using Example 15.4.1 from a free calculus text/15%3A_Multiple_Integration/15.4%3A_TripleIntegrals). This example indicates that $$\int{-1}^5 \int_2^4 \int_0^1 x+yz^2~dz~dy~dx=36.$$
x_integral = CompositeRule(-1, 5, 10, MidpointRule)
y_integral = CompositeRule(2, 4, 10, MidpointRule)
z_integral = CompositeRule(0, 1, 10, MidpointRule)
yz_integral = IteratedIntegral(y_integral, z_integral)
xyz_integral = IteratedIntegral(x_integral, yz_integral)
def f(v):
x,y,z = v
return x + y*z**2
xyz_integral.integrate(f)
Integrals over other shapes¶
Iterated integrals are the secret to integrating over other regions. For example, we can integrate $f(x,y)$ over the unit disk $D$ using polar coordinates: $$\iint_D f~dA = \int_0^1 \int_0^{2 \pi} f(r \cos \theta, r \sin \theta)~r~d \theta ~dr.$$ Here is a class which does this:
class DiskIntegral(IntegrationMethod):
def __init__(self, BasicIntegrationMethod, n):
"""Construct the integral over the unit disk.
A BasicIntegrationMethod should be a function or class that returns
an integral approximate over the interval [a, b] when called with the
endpoints a and b.
We compute integrals by a double integral carried out by subdividing
the intervals [0, 1] and [0, 2*pi] into n intervals and then computing
an interated integral.
"""
outer_integral_method = CompositeRule(0.0, 1.0, n, BasicIntegrationMethod)
inner_integral_method = CompositeRule(0.0, 2*np.pi, n, BasicIntegrationMethod)
integral_method = IteratedIntegral(outer_integral_method,
inner_integral_method)
self._im = integral_method
def domain_dimension(self):
"""Return the dimension of the domain of functions being integrated."""
return 2
def integrate(self, f):
"""Return the approximation of the integral of f."""
# we integrate g(r, theta) = f(r cos(theta), r sin(theta))*r
def g(v):
r, theta = v
return f(np.array([r*np.cos(theta), r*np.sin(theta)])) * r
return self._im.integrate(g)
Integrating the function $f(x,y)=2\sqrt{1-x^2-y^2}$ will produce the volume of the unit sphere in $3$-dimensional space. This is $\frac{4}{3} \pi.$
exact_integral = 4/3*np.pi
exact_integral
def f(v):
x,y = v
# maximum was included below because we were getting numerical errors
# near the boundary of the unit circle that were leading 1-x**2-y**2
# to be negative, which then taking a square root was causing errors.
return 2 * np.sqrt( np.maximum(0.0, 1-x**2-y**2) )
disk_integral = DiskIntegral(MidpointRule, 10)
disk_integral.integrate(f)
disk_integral = DiskIntegral(TrapezoidRule, 10)
disk_integral.integrate(f)
disk_integral = DiskIntegral(SimpsonsRule, 5)
disk_integral.integrate(f)
Remarks:
In the spirit of writing more general code, it might be good to write a class that handles a general change of coordinates in integration. I didn't do this because it would have required computing the Jacobian (derivative matrix).
For a more serious software package for numerically integrating over interesting regions (and with respect to natural measures) see quadpy.
