CCL/CPL

APPENDIX 1

Making Shadow Classes with SWIG

A shadow class in Python is really a Python class. It calls various C++ functions generated by SWIG, which in turn call the original C++ functions. A template class such as mac::vector3D<T> needs to be slightly modified before using SWIG. This appendix shows some examples on how to modify the original C++ code and to instruct SWIG by writing the interface file.

A1.1 Modifying the Source Code

As described above, the first step is to copy the source code in C++ to another place, and make necessary changes. It is still hard to write the C++ code directly compatible with SWIG. These changes include the following.

Nested class

Nested class definition must be removed or commented out before using SWIG, since it is not supported.

Overloading

The order of overloaded functions and operators may be sorted such that the one will appear in Python is the first one. The later ones will be ignored. If they are not sorted, explicit instruction is needed in the interface file in order to ignore some overloaded functions and operators. If overloaded functions and operators are desired in Python, some of them must be renamed in the interface file. See below for more.

Non-member functions

Unbounded non-member functions should be removed. If they are really needed, they should reside in a separate function module. I have not found a way to merge functions with classes in one module.

Additional member functions

Additional member functions specific to Python shadow class may be included, such as __str__, __repr__, etc.

Others

It may be desirable to merge the .c++ file into the header file, but if the .c++ file is too long, it can be included at the end of the header file by using #include, and still kept as a separate file. Other unnecessary bits can be removed in order to make the SWIG interface file and output log cleaner.

The following listing is the modified vector3D.h.

#ifndef CPL_MAC_VECTOR3D_H

#define CPL_MAC_VECTOR3D_H

#include <iostream.h>

#include <strstream.h>

#include <fstream.h>

#include <algo.h>

#include <complex.h>

#include <float.h>

#include "utility.h"

#include "matrix3D.h"

namespace mac

{

template <class T>

class vector3D

{

friend ostream &operator<<(ostream& out, const vector3D<T>& v)

{out << v.getX() << "," << v.getY() << "," << v.getZ();

return out;}

public:

vector3D(const T u = 0, const T v = 0, const T w = 0)

: x(u), y(v), z(w) {}

vector3D(const vector3D<T>& v) : x(v.x), y(v.y), z(v.z) {}

~vector3D() {}

// Conversion

// operator vector3() const {return vector3(x, y, z);}

// operator vector3D<double>() const {return vector3D<double>(x, y, z);}

// operator vector3D<int> () const

// {return vector3D<int>(int(x), int(y), int(z));}

vector3D<int> nearest() const

{return vector3D<int>(mac::nearest(x),

mac::nearest(y),

mac::nearest(z));}

// Output

// void write(ofstream* f) const {if (f) *f << *this << endl;}

// void write(const char* = message::typeGeneral) const;

#ifdef SWIG

%extend

{char *__str__()

{ostrstream line;

line << "<" << self->getX() << " "

<< self->getY() << " "

<< self->getZ() << ">" << ends;

return line.str();}}

#endif

// Element access

inline T operator()(const int i) const

{switch (i) {case 1: return x; case 2: return y; case 3: return z;

default: assert (0);}

return x;} // This will never happen.

inline T& operator()(const int i) // 1-offset, i = 1 - 3

{switch (i) {case 1: return x; case 2: return y; case 3: return z;

default: assert (0);}

return x;} // This will never happen.

inline T operator[](const int i) const {return (*this)(i + 1);}

// 0-offset, i = 0 - 2

inline T& operator[](const int i) {return (*this)(i + 1);}

T getX() const {return x;}

T getY() const {return y;}

T getZ() const {return z;}

// Assignment

const vector3D<T> &operator=(const vector3D<T>& v)

{if (&v != this) {x = v.x; y = v.y; z = v.z;} return *this;}

// Addition/subtraction

vector3D<T> operator+=(const vector3D<T>& v)

{x += v.x; y += v.y; z += v.z; return *this;}

vector3D<T> operator-=(const vector3D<T>& v)

{x -= v.x; y -= v.y; z -= v.z; return *this;}

vector3D<T> operator+(const vector3D<T>& v) const

{vector3D<T> u = *this; return u += v;}

vector3D<T> operator-(const vector3D<T>& v) const

{vector3D<T> u = *this; return u -= v;}

vector3D<T> operator+() const {return *this;}

vector3D<T> operator-() const {return vector3D<T>(0, 0, 0) - *this;}

// Pre-multiplication

vector3D<T> operator*=(const mac::matrix3D<T>& m)

{vector3D<T> v = *this;

x = m(1, 1) * v.x + m(1, 2) * v.y + m(1, 3) * v.z;

y = m(2, 1) * v.x + m(2, 2) * v.y + m(2, 3) * v.z;

z = m(3, 1) * v.x + m(3, 2) * v.y + m(3, 3) * v.z;

return *this;}

vector3D<T> operator/=(const mac::matrix3D<T>& m)

{mac::matrix3D<T> n = m; n.inverse(); return *this *= n;}

vector3D<T> operator*(const mac::matrix3D<T>& m) const

// Caution: This operator means matrix3D<T> * vector3D<T>, returns a

// vector3D<T>.

{vector3D<T> v = *this; return v *= m;}

// Scaling

vector3D<T> operator*=(const T a) // isotropic scaling

{x *= a; y *= a; z *= a; return *this;}

vector3D<T> operator/=(const T a) // isotropic scaling

{x /= a; y /= a; z /= a; return *this;}

vector3D<T> operator*=(const vector3D<T>& v)

{x *= v.x; y *= v.y; z *= v.z; return *this;} // anisotropic scaling

vector3D<T> operator/=(const vector3D<T>& v)

{assert (v.x && v.y && v.z);

x /= v.x; y /= v.y; z /= v.z;

return *this;} // anisotropic scaling

vector3D<T> operator*(const T a) const // isotropic scaling

{vector3D<T> v = *this; return v *= a;}

vector3D<T> operator/(const T a) const // isotropic scaling

{vector3D<T> v = *this; return v /= a;}

vector3D<T> operator/(const vector3D<T>& v) const // anisotropic scaling

{vector3D<T> u = *this; return u /= v;}

// Cross (vector) multiplication

static vector3D<T> cross(const vector3D<T>& u, const vector3D<T>& v)

{return vector3D<T>(u.y * v.z - u.z * v.y,

u.z * v.x - u.x * v.z,

u.x * v.y - u.y * v.x);}

vector3D<T> cross(const vector3D<T> &v) const {return cross(*this, v);}

vector3D<T> operator->*(const vector3D<T> &v) const {return cross(v);}

// Dot (scalar) multiplication

static T dot(const vector3D<T> &u, const vector3D<T> &v)

{return u.x * v.x + u.y * v.y + u.z * v.z;}

T dot(const vector3D<T> &v) const {return dot(*this, v);}

T operator*(const vector3D<T> &v) const {return dot(v);}

// Length

static double vabs(const T x, const T y, const T z)

{return sqrt(double(x * x + y * y + z * z));}

static double vabs(const vector3D<T> &v) {return vabs(v.x, v.y, v.z);}

double vabs() const {return vabs(*this);}

static T vnorm(const T x, const T y, const T z)

{return vnorm(vector3D<T>(x, y, z));}

static T vnorm(const vector3D<T> &v)

{return v.x * v.x + v.y * v.y + v.z * v.z;}

T vnorm() const {return vnorm(*this);}

static bool shorterThan(const vector3D<T>& u, const vector3D<T>& v)

{const T un = u.vnorm();

const T vn = v.vnorm();

if (un == vn) return u < v;

else return un < vn;}

bool shorterThan(const vector3D<T> &v) const

{return shorterThan(*this, v);}

// class shorter

// {

// public:

// bool operator()(const vector3D<T>& u, const vector3D<T>& v) const

// {return shorterThan(u, v);}

// };

// Comparison

bool operator!() const {return *this == vector3D<T>(0, 0, 0);}

bool operator==(const vector3D<T> &v) const

{return x == v.getX() && y == v.getY() && z == v.getZ();}

bool operator!=(const vector3D<T> &v) const

{return x != v.getX() || y != v.getY() || z != v.getZ();}

bool operator&=(const vector3D<T> &v) const // approximately equal

{return vnorm(*this - v) < 10 * FLT_EPSILON;}

bool operator> (const vector3D<T> &v) const

{return !(*this < v) && !(*this == v);}

bool operator< (const vector3D<T> &v) const

{return lessThan(*this, v);}

bool operator>=(const vector3D<T> &v) const

{return *this > v || *this == v;}

bool operator<=(const vector3D<T> &v) const

{return *this < v || *this == v;}

static bool lessThan(const vector3D<T>& u, const vector3D<T>& v)

{if (u.x == v.x)

if (u.y == v.y) return u.z < v.z;

else return u.y < v.y;

else return u.x < v.x;}

// class less

// {

// public:

// bool operator()(const vector3D<T>& u, const vector3D<T>& v) const

// {return lessThan(u, v);}

// };

// Projection

double projectionOn(const vector3D<T>& v) const

{assert (!!v); return dot(v) / v.vabs();}

// Projection of *this on v.

complex<double> projectionOn(const vector3D<T>& u, const vector3D<T>& v)

const {return complex<double>(projectionOn(u), projectionOn(v));}

// Projection of *this on u and v.

vector3D<double> projectionOn(const vector3D<T>& u,

const vector3D<T>& v,

const vector3D<T>& w) const

{return vector3D<double>(projectionOn(u),

projectionOn(v),

projectionOn(w));}

// Projection of *this on u, v, and w.

double projectionOf(const vector3D<T>& v) const

{assert (!!*this); return dot(v) / vabs();}

// Projection of v on *this.

static double projection(const vector3D<T>& u, const vector3D<T>& v)

{assert (!!v); return u * v / v.vabs();} // Projection of u on v.

static T projectionSq(const vector3D<T>& u, const vector3D<T>& v)

{assert (!!v);

const T d = u * v;

const T p = d * d / v.vnorm();

return d >= 0 ? p : -p;}

// The squared projection of u on v. Do not use with T being any

// integer type.

// Orientation

vector3D<double> vcos() const

{assert (!!*this); const double r = vabs();

return r ? vector3D<double>(x, y, z) / r : vector3D<double>(0, 0, 0);}

//return r ? *this / r : vector3D<double>(0, 0, 0);}

vector3D<double> vsin() const

{const vector3D<double> a = orientationAngle();

return vector3D<double>(sin(a.getX()),

sin(a.getY()),

sin(a.getZ()));}

vector3D<double> vtan() const

{const vector3D<double> a = orientationAngle();

return vector3D<double>(tan(a.getX()),

tan(a.getY()),

tan(a.getZ()));}

vector3D<double> orientationAngle() const

{const vector3D<double> c = vcos();

return vector3D<double>(acos(c.getX()),

acos(c.getY()),

acos(c.getZ()));}

complex<double> polarAngle() const

// Polar angles psi and phi define an orientation. In a orthogonal

// coordinate system, rotating Y-axis around Z-axis by psi according

// to the left-hand rule, followed by a rotation of the new Y-axis

// around the old Y-axis by phi according to the right-hand rule, the

// resulting Y-axis is the orientation defined by polar angles psi

// and phi (Z. Ren, Thesis, Academia Sinica, Beijing, China, 1990).

{const vector3D<double> c = vcos();

const double psi = acos(c.getY());

// 0 <= psi <= pi, therefore sin(psi) >= 0.

double phi; // 0 <= phi <= pi, so the return phi is (-pi, pi].

const double sinPsi = sin(psi);

if (!sinPsi) phi = 0;

else {const double cosPhi = c.getX() / sinPsi;

if (cosPhi >= 1) phi = 0;

else if (cosPhi <= -1) phi = M_PI;

else phi = acos(cosPhi);}

return c.getZ() > 0 ? complex<double>(psi, -phi) :

complex<double>(psi, phi);}

vector3D<double> polarAngle(const double kappa) const

{const complex<double> polar = polarAngle();

return vector3D<double>(real(polar), imag(polar), kappa);}

bool operator||(const vector3D<T> &v) const

{return this->vcos() == v.vcos() && !!v;}

bool parallel(const vector3D<T> &v) const {return *this || v;}

bool antiparallel(const vector3D<T> &v) const {return *this || -v;}

bool perpendicular(const vector3D<T> &v) const {return !(*this * v);}

// Distance

static double distance(const vector3D<T> &u, const vector3D<T> &v)

{return vabs(u.x - v.x, u.y - v.y, u.z - v.z);}

double distance(const vector3D<T> &v) const {return distance(*this, v);}

double operator|(const vector3D<T> &v) const {return distance(v);}

// Angle

static double cosAngle(const vector3D<T> &u, const vector3D<T> &v)

{return (u * v) / (vabs(u) * vabs(v));}

double cosAngle(const vector3D<T> &v) const {return cosAngle(*this, v);}

static double sinAngle(const vector3D<T> &u, const vector3D<T> &v)

{return vabs(u ->* v) / (vabs(u) * vabs(v));}

double sinAngle(const vector3D<T> &v) const {return sinAngle(*this, v);}

static double tanAngle(const vector3D<T> &u, const vector3D<T> &v)

{return vabs(u ->* v) / (u * v);}

double tanAngle(const vector3D<T> &v) const {return tanAngle(*this, v);}

static double angle(const vector3D<T> &u, const vector3D<T> &v)

// in radian

{return atan2(sinAngle(u, v), cosAngle(u, v));}

double angle(const vector3D<T> &v) const {return angle(*this, v);}

double operator^(const vector3D<T> &v) const {return angle(v);}

// Area

static double triangleVertices(const vector3D<T> &u,

const vector3D<T> &v,

const vector3D<T> &w)

// u, v, w are vertices

{return vabs((v - u) ->* (w - u)) / 2;}

double triangleVertices(const vector3D<T> &u, const vector3D<T> &v) const

{return triangleVertices(*this, u, v);}

static double triangleEdges(const vector3D<T> &u,

const vector3D<T> &v)

// u and v are two edges

{return triangleVertices(vector3D<T>(0, 0, 0), u, v);}

double triangleEdges(const vector3D<T> &v) const

{return triangleEdges(*this, v);}

// Volume

static T scalarTripleProduct(const vector3D<T> &u,

const vector3D<T> &v,

const vector3D<T> &w)

{return u * (v ->* w);}

T scalarTripleProduct(const vector3D<T> &u, const vector3D<T> &v) const

{return scalarTripleProduct(*this, u, v);}

static T parallelepipedVolume(const vector3D<T> &u,

const vector3D<T> &v,

const vector3D<T> &w)

{const T stp = scalarTripleProduct(u, v, w);

return stp < 0 ? -stp : stp;}

T parallelepipedVolume(const vector3D<T> &u, const vector3D<T> &v) const

{return parallelepipedVolume(*this, u, v);}

// Location

bool inside(const vector3D<T> &u, const vector3D<T> &v) const

// Return whether *this is inside the parallelepiped whose two body

// diagonal vertices are u and v, respectively.

{return min(u.x, v.x) <= x && x <= max(u.x, v.x) &&

min(u.y, v.y) <= y && y <= max(u.y, v.y) &&

min(u.z, v.z) <= z && z <= max(u.z, v.z);}

bool outside(const vector3D<T> &u, const vector3D<T> &v) const

{return !inside(u, v);}

private:

T x, y, z;

};

} // namespace mac

template <class T>

mac::vector3D<double> sin(const mac::vector3D<T>& v) {return v.vsin();}

template <class T>

mac::vector3D<double> cos(const mac::vector3D<T>& v) {return v.vcos();}

template <class T>

mac::vector3D<double> tan(const mac::vector3D<T>& v) {return v.vtan();}

template <class T>

double abs(const mac::vector3D<T>& v) {return v.vabs();}

template <class T>

double abs(const T& x, const T& y, const T& z)

{return mac::vector3D<T>(x, y, z).vabs();}

template <class T>

double norm(const mac::vector3D<T>& v) {return v.vnorm();}

template <class T>

double norm(const T& x, const T& y, const T& z)

{return mac::vector3D<T>(x, y, z).vnorm();}

template <class T>

double vabs(const mac::vector3D<T>& v) {return v.vabs();}

template <class T>

T vnorm(const mac::vector3D<T>& v) {return v.vnorm();}

template <class T>

double vabs(const T& x, const T& y, const T& z)

{return mac::vector3D<T>(x, y, z).vabs();}

template <class T>

T vnorm(const T& x, const T& y, const T& z)

{return mac::vector3D<T>(x, y, z).vnorm();}

template <class T>

T dot_prod(const mac::vector3D<T>& v1,

const mac::vector3D<T>& v2) {return v1 * v2;}

template <class T>

mac::vector3D<T> operator*(const T& x,

const mac::vector3D<T>& v) {return v * x;}

template <class T>

mac::vector3D<T> operator/(const T& x, const mac::vector3D<T>& v)

{return mac::vector3D<T>(x, x, x) / v;}

template <class T>

mac::vector3D<T> operator*(const mac::matrix3D<T>& m,

const mac::vector3D<T>& v)

{return v * m;} // pre-multiplication

#endif

Listing A1.1.0.0.1 The modified vector3D.h ready for wrapping by SWIG.

Several shadow classes can be included in a module. Listing A1.1.0.0.2 shows matrix3D.h. This header file includes the .c++ file at the end.

#ifndef CPL_MAC_MATRIX3D_H

#define CPL_MAC_MATRIX3D_H

#include <iostream.h>

#include <strstream.h>

#include <assert.h>

//#include "matrix.h"

namespace mac

{

template <class T> class vector3D;

template <class T>