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
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;
}
/* 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_
%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"
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.
