Software design: interfaces

Today we are going to talk about software: what to do, what not to do, and why it matters. This falls squarely in the category of "stuff I wish someone had pointed out earlier, that way I wouldn't have had to learn it the hard way." To cut to the chase: this'll mostly be about interfaces, and how to design them to make your code maintainable, usable, and extensible. I'll be working with examples in Python, but the core ideas are more or less language-independent.
I've found that not everyone is entirely convinced that this stuff is important, so I want to spend some time up-front addressing that. If you are already convinced that paying a little more attention to software design is a good thing, great! Go ahead and jump to the section where we look at some code. If not, let me try to address some of the more common objections to this stuff, which hopefully might convince you otherwise.

Objection: I need results fast. I don't have time to make my code nice!

...And that's a totally fair point. I get it: we live in a world where clients need solutions quickly; where deadlines loom and an imperfect solution is better than nothing at all. Why waste time on making sure things are nicely wrapped with a bow on top?
The answer, in a nutshell, has to do with (1) how complex your project is, and (2) how robust it needs to be. If all you're doing is writing "Hello World," then none of this applies and you should ignore everything I'm saying. But what about if you're designing a program to compute the airflow over an airplane? What if you're designing a product that has users in mind, and these users can't be expected to be experts? What if someone would like to adapt your code to handle new tasks after you've written it — what if that person is YOU, two years after you wrote it?
I don't need to belabor the point (said everyone who has ever belabored the point, as I will now do). But the take-home is this. If you're doing something very complicated, and if there are people who would like to use what you are doing in an accessible way, great; that means it's important! And in that case, if you don't pay attention to proper software design, what you're doing is, you're creating a giant maze, and you're plunking your users (and likely also yourself) right down in the middle of it, with no food, no water and certainly no signposts on how to get out. You might make Franz Kafka proud, but that's an issue of literary aesthetics, not useful engineering.
And with that fully belabored, let's get on with it.

Code sample #1, a.k.a. "the wrong way"

Alright. Let's start with the "wrong" way to do things. For this post, we'll use the example of a custom solver for ordinary differential equations (ODEs). In real life, you would use a third-party library, such as Scipy in Python, to integrate ODEs, but this is a convenient example for our purposes.
So in the hypothetical world of this post, you've just written this code here: 
          
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d

t      = np.linspace( 0 , 10 , 1000 )
integration_type = 'forward_euler'

dt     = t[1] - t[0]
X , P  = [] , []
params = [ 10 , 28 , 8./3 ] # sigma , rho , beta
x0     = np.random.uniform( [ 0.5 , 0.5 , 0.5 ] )

x      = np.zeros( [ len(x0) , len(t) ] )
x[:,0] = x0

if integration_type == 'forward_euler':

    for i in range(len(t)-1):
        
        x[0,i+1] = x[0,i] + dt * ( params[0] * ( x[1,i] - x[0,i] ) )
        x[1,i+1] = x[1,i] + dt * ( x[0,i] * ( params[1] - x[2,i] ) - x[1,i] )
        x[2,i+1] = x[2,i] + dt * ( x[0,i] * x[1,i] - params[2] * x[2,i] )

elif integration_type == 'adams_bashforth':

    for i in range(len(t)-1):
        
        f_x = ( params[0] * ( x[1,i] - x[0,i] ) )
        f_y = ( x[0,i] * ( params[1] - x[2,i] ) - x[1,i] )
        f_z = ( x[0,i] * x[1,i] - params[2] * x[2,i] )

        f_xyz = np.array( [ f_x , f_y , f_z ] )

        xyz_pred = x[:,i] + dt * f_xyz

        f_xPred = ( params[0] * ( xyz_pred[1] - xyz_pred[0] ) )
        f_yPred = ( xyz_pred[0] * ( params[1] - xyz_pred[2] ) - xyz_pred[1] )
        f_zPred = ( xyz_pred[0] * xyz_pred[1] - params[2] * xyz_pred[2] )

        f_xyzPred = np.array( [ f_xPred , f_yPred , f_zPred ] )

        x[:,i+1] = xyz_pred + 1.5 * dt * f_xyzPred - 0.5 * dt * f_xyz

elif integration_type == 'rk4':

    for i in range(len(t)-1):

        f_x = ( params[0] * ( x[1,i] - x[0,i] ) )
        f_y = ( x[0,i] * ( params[1] - x[2,i] ) - x[1,i] )
        f_z = ( x[0,i] * x[1,i] - params[2] * x[2,i] )
        k1 = np.array( [ f_x , f_y , f_z ] )

        f_x_k2 = ( params[0] * ( ( x[1,i] + dt/2 * k1[1] ) - ( x[0,i] + dt/2 * k1[0] ) ) )
        f_y_k2 = ( ( x[0,i] + dt/2 * k1[0] ) * ( params[1] - ( x[2,i] + dt/2 * k1[2] ) ) - ( x[1,i] + dt/2 * k1[1] ) )
        f_z_k2 = ( ( x[0,i] + dt/2 * k1[0] ) * ( x[1,i] + dt/2 * k1[1] ) - params[2] * ( x[2,i] + dt/2 * k1[2] ) )
        k2 = np.array( [ f_x_k2 , f_y_k2 , f_z_k2 ] )

        f_x_k3 = ( params[0] * ( ( x[1,i] + dt/2 * k2[1] ) - ( x[0,i] + dt/2 * k2[0] ) ) )
        f_y_k3 = ( ( x[0,i] + dt/2 * k2[0] ) * ( params[1] - ( x[2,i] + dt/2 * k2[2] ) ) - ( x[1,i] + dt/2 * k2[1] ) )
        f_z_k3 = ( ( x[0,i] + dt/2 * k2[0] ) * ( x[1,i] + dt/2 * k2[1] ) - params[2] * ( x[2,i] + dt/2 * k2[2] ) )
        k3 = np.array( [ f_x_k3 , f_y_k3 , f_z_k3 ] )

        f_x_k4 = ( params[0] * ( ( x[1,i] + dt/2 * k3[1] ) - ( x[0,i] + dt/2 * k3[0] ) ) )
        f_y_k4 = ( ( x[0,i] + dt/2 * k3[0] ) * ( params[1] - ( x[2,i] + dt/2 * k3[2] ) ) - ( x[1,i] + dt/2 * k3[1] ) )
        f_z_k4 = ( ( x[0,i] + dt/2 * k3[0] ) * ( x[1,i] + dt/2 * k3[1] ) - params[2] * ( x[2,i] + dt/2 * k3[2] ) )
        k4 = np.array( [ f_x_k4 , f_y_k4 , f_z_k4 ] )

        x[:,i+1] = x[:,i] + 1./6 * dt * ( k1 + 2*k2 + 2*k3 + k4 )
    

fig = plt.figure()
ax  = fig.add_subplot( projection='3d' )
ax.plot3D( x[0] , x[1] , x[2] )
plt.show()
What's wrong with this? In the strict sense of functionality, nothing. It does exactly what you designed it to do: integrate the Lorenz equations, and it even makes a nice plot of the solution.
Excited by your success, your client says, "great! Could you also make it integrate the Duffing equation?"
And that's where the wind leaves your sails.

Refactoring: reducing dependencies

The biggest obstacle to fulfilling the request of your client — who has a strange affinity for toy dynamical systems, it would seem — is that the code which implements the time integrator is completely entangled with the code that implements the right-hand side. This makes it impossible to implement a new right-hand side (Duffing) and "plug it into" your time integrator (e.g., RK4) without either (1) destroying the code you already wrote and replacing it, or (2) duplicating code for the integrator, along with a horror show of even more kludgy "if-else" conditionals to mediate which right-hand side is integrated.
The solution here is twofold: first, we need to separate the right-hand side from the integrator; doing so will make it much easier to implement new right-hand sides while leaving the integrator code unchanged. Of course, these are not completely independent concepts — after all, a time integrator exists to integrate a right-hand side. Recognizing this leads to the second part of the solution: we should link the integrator and right-hand side objects together through composition, which will allow us to implement the loose dependencies we are looking for.
By now, you may be starting to have a high-level vision that looks something like this: 
              
# make the RHS object first, with the parameters needed to define it
rhs        = Lorenz( { "sigma" : 10 , "rho" : 28 , "beta" : 8./3 } )
# integrator is composed of a rhs (and possibly other things)
integrator = Forward_Euler( rhs , [ possibly other args ] )
We're ready to take a look at the full version of this below.

Code sample #2, a.k.a. "the right way"

This is the high-level refactored version of the original ODE solver:
          
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits import mplot3d
from integrator import *
from ode_rhs import *

rhs        = Lorenz( { "sigma" : 10 , "rho" : 28 , "beta" : 8./3 } )
integrator = Forward_Euler( rhs , dt=0.01 , t_steps=1000 )
x          = integrator.solve( np.array( [ 0.5 , 0.5 , 0.5 ] ) )

fig = plt.figure()
ax  = fig.add_subplot( projection='3d' )
ax.plot3D( x[0] , x[1] , x[2] )
plt.show()
One thing that immediately jumps out when comparing this code with the other is that this version is shorter and much easier to understand. Why is that? For two reasons: one, as we've already discussed, we've disentangled the right-hand side and time integrator; they are now distinct pieces of code. But the biggest reason this is easier on the eyes is because we've hidden the implementation details. Remember that this is just the highest level code: nowhere does it say here how Lorenz or Forward_Euler actually work. That's intentional: the user shouldn't have to care about that.

Implementation details: time integrator

There are two interesting questions we need to answer regarding the time integrator: what defines a valid time integrator, and how does it use its right-hand side (given to it at construction)? The first question is important because there are obviously many schemes for time integration (e.g., forward Euler, Runge Kutta). In order for all of them to be usable in consistent ways, we'd like to unite them with inheritance. Many would argue (as would I) that inheritance works best when it's shallow and broad; this prevents all kinds of dangerous things from happening (see: Diamond inheritance).
With that in mind, we make a very basic base class that defines the minimal set of things that a time integrator needs to do, like so: 
          
class Integrator(ABC):

    def __init__( self , 
                  rhs : RHS , 
                  dt : float , 
                  t_steps : int ):
        
        """
        Base class for time integrators.

        Args: 
            rhs: instance of RHS
            dt: float, desired (fixed) time-step
            t_steps: int, desired number of time-steps
        """

        self._rhs  = rhs
        self._dt   = dt
        self._n    = t_steps

    @abstractmethod
    def solve( self , 
                x0 : np.ndarray ) -> np.ndarray:
        
        """
        Integrate the ODE forward in time.
        
        Args:
            x0: np.ndarray, initial state values
        Returns:
            2-D numpy array of size (X,T) where X is state dimension and T is number of time steps
        """

        pass
As we can see, an Integrator is something that owns a RHS object (to be defined), along with a desired timestep and final integration time. It also must implement the method Integrator::solve(x0), which maps a given initial condition x0 to a matrix containing the solution in time. And that's it. Perhaps the take-away of this section is how minimal the requirements are for what an Integrator must do. Again, these are implementation details and so are someone else's problem to worry about.
We'll finish here by giving one concrete example implementation of Integrator (for forward Euler): 
          
class Forward_Euler( Integrator ):

    def solve( self , 
                x0 : np.ndarray ) -> np.ndarray:

        x      = np.zeros( ( len(x0) , self._n ) )
        x[:,0] = x0

        for i in range( self._n-1 ):
            x[:,i+1] = x[:,i] + self._dt * self._rhs( x[:,i] )

        return x

Implementation details: right-hand side

The design of this class hierarchy should be obvious given what we've asked it to be capable of doing in the previous higher-level code snippets. In words, a right-hand side should be an object that is constructed with some parameter name/value pairs, and which a user can call to output the right-hand side evaluation at a given point. Translated into code, we have this: 
              
class RHS(ABC):
    """
    base class for a right-hand side
    """

    def __init__( self , param : Dict[ str , float ] ):
        """
        inputs:
            param : dict mapping parameter names to values
        """
        self._param = param

    @abstractmethod
    def __call__( self , 
                  x : np.ndarray ) -> np.ndarray:

        pass
Here's one example implementation of RHS: 
              
class Lorenz( RHS ):

    def __call__( self , 
                  x : np.ndarray ) -> np.ndarray:
        
        """
        inputs
            x : initial state

        example : 
            >>> params = { "sigma" : 10 , "rho" : 28 , "beta" : 8./3 }
            >>> rhs    = Lorenz( params )
            >>> x0     = np.array( [1,1,1] )
            >>> dxdt   = rhs( x0 )
        """

        f_x = self._param['sigma'] * ( x[1] - x[0] )
        f_y = x[0] * ( self._param['rho'] - x[2] ) - x[1]
        f_z = x[0] * x[1] - self._param['beta'] * x[2]
        
        return np.array( [ f_x , f_y , f_z ] )

Parting thoughts

Hopefully it should be clear now why the second approach is preferable to the initial one: recall our hypothetical client's request of extending the code's functionality to cover a brand new use case (i.e., a new right-hand side). Previously, this was not possible without taking a sledgehammer to everything we'd done, since the code for our time integrator was tightly coupled to the code for the right-hand side. Now, we've largely decoupled those two things. Extending the functionality of either the right-hand side or the time integrator is now a process that is easy to do without one affecting the other, and moreover, there are clear guidelines (through inheritance) on how to do these things.
Another thing to consider: the same decoupling that made this code easier to extend has also made it much better from the standpoint of test coverage. If you were to encounter bugs in running the first version of the code, it might not be immediately clear whether those errors stemmed from an problem with the time integrator, or with the right-hand side (or both). Now, because these are separate objects, it is possible to write unit tests to verify the functionality of both, independently.