SWIG and Python | chriscoughlin.com

I originally wrote this up earlier this year as part of a how-to on creating plugins for NDIToolbox. I thought it might be worth posting it on its own in case anyone’s looking for a quick intro to using SWIG with Python and C++. Among other things, SWIG is great if you started with a pure Python program but find yourself needing a bit of a speed boost: just take the number-crunching or other expensive code and replace it with C or C++.

Pure Python Application

To demonstrate how to build and structure a Python-only toolkit for NDIToolbox we’ll build a Helmholtz coils calculator to calculate the magnetic field produced by a pair of identical coils of radius a, number of turns N, and each carrying a current I. If the coils are in opposition (either the coils are wound in opposite directions or the coils’ currents are in opposite direction), their magnetic fields interact to create a maximum at the center.

Organization

Following our recommended approach, we’re going to develop our main logic and user interface separately from the toolkit interface. We’ll have one Python file for the logic hhcoils.py that has a class HelmholtzCoils in which all our calculations are contained. We’ve decided to do a simple wxPython user interface, hhcoils_ui.py, that’s responsible for taking user input and reporting the results of our calculator back to the user. Finally, our toolkit interface with NDIToolbox will be in a third Python file hhcoils_toolkit.py.

(Note-I’m skipping the toolkit for this post, and the UI’s been replaced with a simple script that demonstrates function.)

Logic

Our hhcoils.py file that handles all the calculations is as follows. By itself it’s not very interesting but it helps to illustrate a simple way of structuring your application as both a standalone and as a toolkit.

#!/usr/bin/env python

''' HelmholtzCoils - simple magnetic field calculator (originally used to demo plugin development in NDIToolbox)

Chris Coughlin (TRI/Austin, Inc.)
'''

import math

class HelmholtzCoils(object):
    '''Calculates the magnetic field produced by a pair of Helmholtz coils wound in opposition'''
    mu_0 = 4.*math.pi*1e-7

    def __init__(self, turns_per_coil, current_per_coil, coil_radius):
        self.N = int(turns_per_coil) # Turns per coil
        self.I = float(current_per_coil) # Current (A) per coil
        self.a = float(coil_radius) # Common radius (m) of coils
        self.lhcoil_position = -self.a/2. # Position of left coil (m), defined as a/2
        self.rhcoil_position = self.a/2. # Position of right coil (m), defined as a/2
    
    def geometry_correction(self, coil_position, position):
        '''Geometry correction for magnetic field calcs'''
        return math.pow(1 + math.pow(coil_position - position, 2) / math.pow(self.a, 2), -1.5)

    def H(self, position):
        '''Magnetic field at position (m) in A/m'''
        lh_geometry_factor = self.geometry_correction(self.lhcoil_position, position)
        rh_geometry_factor = self.geometry_correction(self.rhcoil_position, position)
        total_geometry_factor = lh_geometry_factor + rh_geometry_factor
        return ((self.N * self.I)/(2 * self.a)) * total_geometry_factor
        
    def centerH(self):
        '''Magnetic field at dead center of coils'''
        return self.H(0)
        
    def B(self, position):
        '''Flux density at position (m) in mT'''
        return HelmholtzCoils.mu_0 * self.H(position)*1000
        
    def centerB(self):
        '''Flux density at dead center of coils'''
        return self.B(0)
        
    def B_mG(self, position):
        '''Flux density at position (m) in mG'''
        return self.B(position) * 1e4
        
    def wirelength(self):
        '''Length of wire (m) to make N turns of radius a'''
        coil_circumference = math.pi * 2. * self.a
        wire_per_coil = float(self.N) * coil_circumference
        return 2.*wire_per_coil
        
    def awg_recommendation(self):
        '''Returns a recommendation for AWG wire gauge to use for the coils for currents between
        0.0125 and 15 Amps.  Conservative estimate, returns 0 as recommendation for currents 
        outside this range.'''
        gauge_current = {10:15, 11:10, 14:5, 17:2.5, 20:1.5, 21:1.0, 22:0.75, 24:0.5,
            27:0.25, 30:0.125, 31:0.100, 32:0.075, 34:0.050, 37:0.025, 40:0.0125}
        delta = float('inf')
        epsilon = 0.
        awg_rec = 0
        if self.I >= 0.0125 and self.I <= 15:
            for gauge, current in gauge_current.items():
                epsilon = current - self.I
                if epsilon >= 0 and epsilon < delta:
                    delta = epsilon
                    awg_rec = gauge
        return awg_rec

            view raw
            hhcoils.py
            This Gist brought to you by GitHub.
          

To use the calculator in our user interface, we’ll create an HelmholtzCoils instance with the user’s choice for the number of turns N per coil, the current per coil I, and the coil radius (and the coil separation distance) a. So by way of a demonstration, let’s pretend the following is your UI…

import hhcoils

if __name__ == "__main__":
    coils = hhcoils.HelmholtzCoils(100, 0.6, 0.8382)
    print("Field at centre (A/m) = {0}".format(coils.centerH()))
    print("Recommend {0} metres of {1} AWG wire.".format(coils.wirelength(),
        coils.awg_recommendation()))

Switching To C++

As an example of how to use more than one language in developing your toolkit, suppose we’ve decided to swap out our Python Helmholtz coils engine for a new high-performance C++ engine. We’re happy with the user interface as-is and the toolkit interface works well; since we kept our layers separate we can simply swap the backends without any other changes.

Logic

We’ve decided that we need a high-performance C++ backend to increase the speed of our calculations. Investigating our options for interfacing C++ and Python we’ve decided to use SWIG. To use SWIG, we’ll first write our C++ Helmholtz coils calculator code; once we’re satisfied with it we’ll write a simple configuration file and let SWIG create a Python wrapper for us.

We’ve decided to have a single simple class HelmholtzCoils, defined in a header hhcoils.h and an implementation file hhcoils.cxx. Notice that we’re using exactly the same method names and arguments as our original Python engine – by doing this the changes are transparent to the GUI and we won’t have to make changes to hhcoils_ui.py. Since our GUI won’t be altered we also won’t have to update our toolkit interface, meaning that we can just drop the new engine into the hhcoils folder and the upgrade is complete. Current customers of our toolkit can do the same with their installations to take advantage of the new engine, or we could write a simple installer to automate the process for them.

/* HelmholtzCoils - demonstrating separation of analysis and toolkit interface in NDIToolbox by creating a
simple magnetic field calculator

Chris Coughlin (TRI/Austin, Inc.)
*/
#include "hhcoils.h"

const double HelmholtzCoils::H(double position) const{
    double lh_geometry_factor = geometry_correction(lhcoil_position,position);
    double rh_geometry_factor = geometry_correction(rhcoil_position,position);
    double total_geometry_factor = lh_geometry_factor + rh_geometry_factor;
    return ((N*I)/(2*a))*total_geometry_factor;
}

const double HelmholtzCoils::geometry_correction(const double coil_position, const double position) const{
    return pow(1 + pow(coil_position - position,2) / pow(a,2), -1.5);
}

const double HelmholtzCoils::wirelength(void) const {
    double coil_circumference = M_PI*2*a;
    double wire_per_coil = N*coil_circumference;
    return 2*wire_per_coil;
}

const double HelmholtzCoils::awg_recommendation(void) const {
    double wire_gauges[] = {10,11,14,17,20,21,22,24,27,30,31,32,34,37,40};
    double currentrecs_awg[] = {15,10,5,2.5,1.5,1.0,0.75,0.5,0.25,0.125,0.100,0.075,0.050,0.025,0.0125};
    double delta = DBL_MAX;
    double epsilon = 0;
    double awg_rec = 0;
    if(I>=0.0125 && I<=15) {
        for(int iter=0;iter<16;iter++) {
            epsilon = currentrecs_awg[iter] - I;
            if(epsilon>=0 && epsilon<delta) {
                delta = epsilon;
                awg_rec = wire_gauges[iter];
            }
        }
    }
    return awg_rec;
}

            view raw
            hhcoils.cxx
            This Gist brought to you by GitHub.
          

/* HelmholtzCoils - demonstrating separation of analysis and toolkit interface in NDIToolbox by creating a
simple magnetic field calculator

Chris Coughlin (TRI/Austin, Inc.)
*/

#ifndef HELMHOLTZCOILS_H_
#define HELMHOLTZCOILS_H_

#define _USE_MATH_DEFINES // Required on some platforms to get mathematical constants
#include <cmath>
#include <cfloat>

static const double mu_0 = 4*M_PI*1e-7;

class HelmholtzCoils {
public:
    HelmholtzCoils(int turns_per_coil, double current_per_coil, double coil_radius):
        N(turns_per_coil), I(current_per_coil), a(coil_radius), lhcoil_position(-a/2), rhcoil_position(a/2) { }
    
    const double H(const double position) const; // Magnetic field at position (m) in A/m
    const double centerH(void) const { return H(0); } // Magnetic field at dead center of coils
    const double B(const double position) const { return mu_0*H(position)*1000; } // Flux density at position (m) in mT
    const double centerB(void) const { return B(0); } // Flux density at dead center of coils
    const double B_mG(const double position) const { return B(position) * 1e4; } // Flux density at position (m) in mG
    
    const double wirelength(void) const; // Length of wire (m) to make N turns of radius a
    const double awg_recommendation(void) const; // Lookup table to make AWG recommendations based on current I
        
private:
    int N; // Turns per coil
    double I; // Current (A) per coil
    double a; // Common radius (m) of coils
    double lhcoil_position; // Position of left coil (m), defined as -a/2
    double rhcoil_position; // Position of right coil (m), defined as a/2
    // Geometry correction for magnetic field calcs
    const double geometry_correction(const double coil_position, const double position) const; 
};

#endif // HELMHOLTZCOILS_H_

            view raw
            hhcoils.h
            This Gist brought to you by GitHub.
          

%module hhcoils
%{
/* Includes the header in the wrapper code */
#include "hhcoils.h"
%}
/* Include various STL interfaces - not really needed here but included as a demo*/
%include "std_string.i"
%include "std_vector.i"
/* For vectors it's necessary to specify what type of vector(s) you'll be using-here we're
using the Python "vector_float" as an alias for vector<float>*/
namespace std {
   %template(vector_float) vector<float>;
   %template(vector_vector_float) vector<std::vector<float> >;
};

/* Parse the header file to generate wrappers */
%include "hhcoils.h" 

view raw hhcoils.i This Gist brought to you by GitHub.

SWIG’s role in our engine swap is to compile the C++ code into a library, and then to create new Python code as a wrapper. Our engine code is simple enough that we don’t require much in the way of configuration, but the config file adds handlers for strings and vectors just to illustrate their use.

To create the library depends on your C++ compiler and your operating system, what version of Python you’re using, etc. but essentially it entails running SWIG on your configuration file, compiling your C++ code and the wrapper code SWIG generates, and finally linking the compiled code into a static library (.so on Linux, .lib on Windows, etc.). For example, running SWIG with GCC on Linux on our toolkit engine:

swig -c++ -python hhcoils.i
c++ -fpic -c hhcoils.cxx hhcoils_wrap.cxx -I /usr/include/python2.7/ -lstdc++
c++ -shared hhcoils_wrap.o hhcoils.o -o _hhcoils.so

This creates a new Python wrapper hhcoils.py; this file and the _hhcoils.so static library is the new C++ Helmholtz coils calculation engine. We can then copy these two files to our pre-existing structure (overwriting our previous Python-only hhcoils.py engine), and the engine upgrade is complete. Our toolkit interface is unchanged and loads our GUI; our GUI already points to hhcoils.py and since we kept the same method names and signatures we don’t have to change the GUI either:

#!/usr/bin/env python
import hhcoils

if __name__ == "__main__":
    coils = hhcoils.HelmholtzCoils(100, 0.6, 0.8382)
    print("Field at centre (A/m) = {0}".format(coils.centerH()))
    print("Recommend {0} metres of {1} AWG wire.".format(coils.wirelength(),
        coils.awg_recommendation()))

Wrapping Up

The nice thing about using this approach to develop our application is we can get the best of both worlds. Python is generally a more productive programming language than C++, which means that we get our product out to the end user faster by developing in Python. If we find later on that we need a little speed boost, we can swap out the slow bits with a little C++ (or another language of your choosing) with little or no trouble.

There are many other ways you can speed up your Python application-the Cython project takes a lot of the heavy lifting out of writing C extensions; PyPy uses a JIT compiler to speed up your code, and so on. Generally speaking if you’re worried about performance you should first make sure your Python code is as efficient as you can make it before you start thinking about these other options, but it’s nice to know they’re there if you need ‘em.