We will use a very simple example, quadrature, to illustrate how to program your own classes in SincLib.  This example uses code from SincLib's Quadrature class.

Theory

We reproduce the following theorem from Lund and Bower's book to illustrate the theoretical underpinning of the sinc quadrature method:


Note that we have exponential convergence and a relatively simple formula for computing quadrature.

Implementation

Here's a snippet of the implementation of sinc quadrature:

import scipy
from sinc.base.SincBaseLib import *
    
log = logging.getLogger('SincMethods.Quadrature')
    
class Quadrature( SincBase ):
  
  @staticmethod
  def defaults():
    """ The defaults of the quadrature class. """
    opts = SincBase.defaults() 
    opts['integrand'] = ['i','integrand','The function to integrate',None,'func']
    opts['problem_file'] = ['p','problemFile','The file describing the quadrature to' + \
                            ' perform.',None]
    opts['answer'] = ['w','answer','The true answer to the integral.',0.0, 'expr']
    return opts
        
  def update( self ): 
    """ Updates the integrand-independent variables. """
    super( Quadrature, self ).update()
    self.nodes = cSincLib.genSincNodes( self.psi, self.M, self.N, self.h )
    self.df2 = self.f2( self.nodes )
    
  def compute(self):
    """ Computes the integral of the given integrand using sinc quadrature. """
    integral = self.integrand( self.nodes ) * self.df2
    self.soln = scipy.sum( integral ) * self.h
    return self.soln

  # The rest of the entry is truncated

In order to take advantage of the underlying library, we declare the Quadrature object to be a subclass of the SincBase object.  SincBase picks out the appropriate transforms and sets up the constants appropriate for the problem.

Next, examine the defaults method.  This function defines all the command line parameters and defaults for the class.  The options which should be returned are a dictionary of lists.  The key of the dictionary is the variable name the option will be saved to in the class.  The value list consists of 5 elements: the short name, the long name, the description, default, and variable type.  The first three have to do with command line options.

Notice that we start off by copying all the options from SincBase.  This is often a good idea, because SincBase defines most of the standard constants needed for sinc numerical analysis.

The two other important functions are update and compute.  Update is run when the problem is first created.  It only needs to be run when one of the problem parameters, such as M or the conformal transform, is changed.  Compute needs to be re-run whenever the integrand function changes.

Update defines two quantities.  First, the sinc nodes are generated using the function cSincLib.genSincNodes and the values of $f_2$ at sinc nodes are pre-computed.  There are many useful sinc-related functions defined in sinc.base.sinclib, and also implemented in C for speed in sinc.base.cSincLib.  The function f₂ is defined as in the table here.

Finally, the sinc quadrature formula is:


This value is computed in the function compute.  First, we evaluate the f at the sinc nodes, then multiply it by the values of f2 at the respective sinc nodes; this produces a vector of size (M+N+1):

integral = self.integrand( self.nodes ) * self.df2

Then, we sum up with the function sum, and multiply by h:

self.soln = scipy.sum( integral ) * self.h

Discussion

We have broken a "rule" of numerical analysis: we keep the intermediate values in the array integral, instead of "summing as we go".  This is done for multiple reasons:

• (M+N+1) is usually small (<400), resulting in an array of size 3.1 KB.  This value can stay in the cache of even a decade-old 486 CPU!  Hence, we don't have to worry about hitting main memory in this computation.
• Python looping constructs (implemented in python byte code) are extremely slow compared to array function implementation (implemented in C); there is a O(100) difference in speed.
• Array manipulation will allow us to parallelize the computation later with minimal effort due to the structure of the PETSc library.
  
• Even for a function g implemented in C, the cost of calling g once for each element of an array is very high compared to functions which must be called once and can manipulate the whole array.  For example, calling sin(X) for the first 40 elements of an array might be a single processor instruction, while something like this:
  for( int i = 0; i < 40; i++ ) { Y[i] = sin(X[i]); }
  might translate into many function calls.

The rule of thumb for SincLib is:

• Maximize the ratio of array operations compared to python operations.  Never, ever, loop through the array in python!  If a large speedup might be gained by manipulating the array directly, investigate in alternates like SciPy's weave module or implementations in C.
• Wherever possible, use numerical library functions as opposed python functions.  If using sum, make sure you use scipy's sum, not Python's (there is a 10x difference).  Remember the writers of LAPACK - the underlying math library - have spent much more time working on optimizations than you have available, and they're experts in the subject.

High Performance Python

For a thorough discussion of high performance code in Python, refer to http://www.scipy.org/PerformancePython.  The writer of that tutorial considers several methods for solving Laplace's equation in Python.  The pure python implementation has a runtime of 1500 seconds; the implementation using Numeric arrays takes 29.3 seconds.  Embedding C into the python code has a runtime of 2.3 seconds.  Finally, writing a pure C++ version results in a runtime of 2.16 seconds. 

Although the difference between Numeric and embedded C is a factor of about 10, the Numeric implementation takes significantly less time to code.  If you are writing a large library (such as SincLib), it's often advantageous to spend your time coding more algorithms, profile your code, and re-implement the core portions in C.  We believe this provides the best balance of maximizing speed and minimizing coding time.