Monthly Archives: September 2012

C++ debugging with GDB and Valgrind quickstart

As I’ve mentioned before, I’m coming into C++ from Python and Matlab. I’ve programmed C++ before, but it’s been a while, and handling the memory management myself and tracking down seg faults again is taking some getting used to. Here are two things that really helped me out.

Program – GDB
As the scale of the program I’m working on in C++ keeps increasing, debugging it becomes more and more of a pain in the ass to do using the old “add print statements” approach to pinpointing the source of the problem. After one seg fault too many I broke down and started looking at how to use GDB. Turns out, it’s very straightforward to get some basic (and very useful) functionality out of it.

Here are the steps:

  1. Add -g to your compile command, ie g++ -g test.cpp -o test
  2. Type in gdb test to your console
  3. Type run into the gdb console, then when your program breaks, type in backtrace or bt to get the line that your program broke down on, and the trace of function calls leading up to it!

There’s a bunch of other functionality built in to gdb, but I was just looking for a very simple, low overhead introduction to it and this was all I needed to get going. Having that line number can save a ton of time.

Program – Valgrind
Once my program got up and going though, I noticed that for longer simulation runs I got core memory dump errors. After a very short search I came across a program called Valgrind that profiles C++ code to find memory leaks. Turns out I had quite a few.
This program is also extremely easy to use. To install on linux is, of course, just sudo apt-get install valgrind. To use Valgrind to check for memory leaks in your program and get information about where they originated, from the console type:

valgrind --leak-check=full --show-reachable=yes --track-origins=yes ./test

If all is going well, then you should see something like this pop up at the end of the printout:

==9422== HEAP SUMMARY:
==9422==     in use at exit: 0 bytes in 0 blocks
==9422==   total heap usage: 31,127 allocs, 31,127 frees, 858,710 bytes allocated
==9422==
==9422== All heap blocks were freed -- no leaks are possible

It’s beautiful.
Valgrind tracks your allocations and freeing of memory, so if it finds that there are unmatched allocations, or you do something fairly tricky with your handling of memory, you’ll get a printout that lets you know there’s unfreed memory and will give you a trace of function calls leading to the memory that isn’t freed. It’s a great program, for more information they have a quickstart guide here: http://valgrind.org/docs/manual/quick-start.html#quick-start.intro.

As I said before, I’m getting used to memory management in C++ again; here are a couple of things that really perplexed me while I was trying to track down the memory leaks:

Perplexment 1 – Copy constructor on dereferenced pointers
I had been passing around pointers into class constructor, and then assigning them to local variables on the stack (so I thought). The idea was then when the class then was destroyed, the stack had taken over control of the memory and would free automatically. It looked like this:

class classA {
public:
    int a;
    VectorXf b;

    classA(); // constructor
};

classA::classA(int a, VectorXf *b) {
    this->a = a;

    if (b == NULL)
        b = new(VectorXf)(VectorXf::Ones(a));
    this->b = *b;
}

WRONG. Apparently this won’t even work at all most of the time, but the Eigen library has an amazing set of copy constructors that let this slide. What’s happening when you dereference this pointer is that a copy of the VectorXf that b points to is created. Then, the instance you created to pass into this constructor is left unfreed floating around in the ether. Terrible. The fix, is, of course, to make b a pointer to a vector instead. i.e.

class classA {
public:
    int a;
    VectorXf* b;

    classA(); // constructor
};

classA::classA(int a, VectorXf *b) {
    this->a = a;

    if (b == NULL)
        b = new(VectorXf)(VectorXf::Ones(a));
    this->b = b;
}

Now this also holds when you return pointers! If you have some function

VectorXf* classA::func1 getVec();

that is returning a pointer to you, you be damned sure you store that result in a pointer. None of this:

VectorXf x = *func1();

or any of its like or ilk, that is bad news. Also this:

VectorXf x = func1()->reverse();

no no no no no no. Memory leak. Same as dereferencing, leaves you with VectorXfs floating around in the wild.

And now you know. And knowledge is power.

Perplexment 2 – Putting virtual in front of base class deconstructor
The other thing that caught me up was that if you have a base class and inheriting classes set up, the deconstructor of the base class must be declared with virtual, or it isn’t called by the inheriting classes. i.e.

virtual ~baseClass();

Additionally, if you don’t have that virtual line in front of your base class deconstructor, none of the inheriting class deconstructors get called either! This was another one that tripped me up.

Perplexment 3 – delete a, b, …
This does not work! It compiles fine, but having a string of variables to delete in a row does not actually delete the latter variables. Perplexing!

Tagged , ,

Getting function run time in C++ on Linux

I wanted to time how long it took for one of my functions to run yesterday and it ended up taking me a little bit of time to track down how to do this in C++ on a Linux box, so I thought I’d throw up the code here.

#include <sys/time.h>
#include <iostream>
using namespace std; 

int main() { 

    struct timeval start_time; 
    struct timeval end_time; 
    
    gettimeofday(&start_time, NULL);
    aFunction(someParams); // function to be timed
    gettimeofday(&end_time, NULL); 

    // get difference, multiply by 1E-6 to convert to seconds
    float duration = (end_t.tv_sec - start_t.tv_sec) + (end_time.tv_usec - start_time.tv_usec) * 1E-6; 
    
    cout << "duration: " << duration << "s" << endl;
}

And that’s it! Not too much to it, but was a little difficult to track down.

Tagged , ,

Cython Journey Part 2: Including Eigen

As I mentioned in my last Cython Journey post, the reason I’m using Cython is to take my speedy C++ code that does a bunch of matrix operations using the Eigen library, and to make them accessible from Python. Now, I have no desire to make an interface for the Eigen library, all I need is to be able to access the MatrixXf objects from the Eigen library, and convert them into numpy.ndarray objects for manipulation in Python. The logical first step towards this end would then seem to be to successfully compile C++ code that links to the Eigen library.

HOWEVER. After messing around with a bunch of linking problems and the like, I came across code by charanpald that provides an interface for Eigen sparse matrices and sparse matrix operations to Python (https://github.com/charanpald/SpPy). In this code, he links to the Eigen library through a buffer class that extends the Eigen::SparseMatrix class and adds some more functionality. There is very possibly a better way to do this, but I also needed to define some extra functions operating on the Eigen matrices, so defining an extension to the MatrixXf class worked out well. For this example, though, I’ll just use a barebones class extension to keep the size manageable. Here is my class:
cpp_matrixxfpy.h

#ifndef MATRIXXFPY_H
#define MATRIXXFPY_H
#include <Eigen/Dense>
using namespace Eigen;

class MatrixXfPy : public MatrixXf { 
    public: 
        MatrixXfPy() : MatrixXf() { }
        MatrixXfPy(int rows,int cols) : MatrixXf(rows,cols) { }
        MatrixXfPy(const MatrixXf other) : MatrixXf(other) { } 
};
#endif

Like, I said, barebones.
With this now we can go to our .pyx file and create a Cython handle into the MatrixXf object.
matrixxfpy.pyx

cdef extern from "cpp_matrixxfpy.h":
    cdef cppclass MatrixXfPy:
        MatrixXfPy()
        MatrixXfPy(int d1, int d2)
        MatrixXfPy(MatrixXfPy other)
        int rows()
        int cols()
        float coeff(int, int)

Same as for the cpp_test code that we wanted to interface with Python; first thing we do in the .pyx file is redefine the class and variables / functions that we want access to in Python, along with the constructors. In this simple example we’re going to only access three functions from MatrixXf, rows(), cols(), and coeff(int, int). These functions return the number of rows, the number of columns, and the matrix coefficient value at a given index.

So now that we have a Cython handle on the Eigen code, let’s make this functionality available to Python. HOWEVER. Let’s keep in mind the goal, we’re not here to provide a Python wrapper for MatrixXf objects. What we really want is just to have Eigen and our C++ code over there, somewhere else, doing all of the calculations etc and then for us to be over in Python and just be able to get a numpy.ndarray for Python to play with and plot. Let us set about that now.

Let’s go back to the Test class we defined in Cython Journey Part 1. First thing is first, we need to edit the C++ code to create and return a MatrixXf object. So let’s just throw in a simple function:
cpp_test.h

#ifndef TEST_H
#define TEST_H

#include "cpp_matrixxfpy.h"
using namespace Eigen; 

class Test { 
public: 
    int test1; 
    Test();
    Test(int test1); 
    ~Test(); 
    int returnFour(); 
    int returnFive();
    Test operator+(const Test& other); 
    Test operator-(const Test& other);
    MatrixXfPy getMatrixXf(int d1, int d2);    
};
#endif

and define it in cpp_test.cpp as

#include "cpp_test.h"

Test::Test() { 
    test1 = 0;
}

Test::Test(int test1) { 
    this->test1 = test1; 
}

Test::~Test() { }

int Test::returnFour() { return 4; }

int Test::returnFive() { return returnFour() + 1; }

Test Test::operator+(const Test& other) { 
    return Test(test1 += other.test1);
}

Test Test::operator-(const Test& other) { 
    return Test(test1 -= other.test1);
}

MatrixXfPy Test::getMatrixXf(int d1, int d2) { 
    MatrixXfPy matrix = (MatrixXf)MatrixXfPy::Ones(d1,d2);
    matrix(0,0) = -10.0101003; // some manipulation, to show it carries over
    return matrix;
}

I put up the whole class again so it’s all in one place.

Alright, that’s done. Now we can look at our .pyx file for Test. We’re going to need to make the MatrixXfPy class defined just above available to our Test class. To do this is easy, all we need to throw in is:
test.pyx

 
include "cpp_matrixxfpy.h"

at the top! And while we’re at it, we’re going to be adding some more functionality to test that requires numpy, so also throw in:

 
import numpy

This additional functionality is going to access the getMatrix(int,int), create an numpy.ndarray object, and fill it out with the values from the matrix. Let’s get about doing that. We need to expose the new C++ function to Cython, so add this line to the cdef cppclass Test declarations:

        MatrixXfPy getMatrixXf(int d1, int d2)

And now we can add the code that will be accessible from Python. I’ll post the code and then describe it below (although it’s pretty straightforward):

    def getNDArray(self, int d1, int d2): 
        cdef MatrixXfPy me = self.thisptr.returnMatrixXf(d1,d2) # get MatrixXfPy object
        
        result = numpy.zeros((me.rows(),me.cols())) # create nd array 
        # Fill out the nd array with MatrixXf elements 
        for row in range(me.rows()): 
            for col in range(me.cols()): 
                result[row, col] = me.coeff(row, col)   
        return result 

I just flat out stripped the form for this from the SpPy project and simplified it, so it’s not going to be optimized for speed in any way, but it gets the job done and it’s very easy to read. First we just set up a numpy.zeros matrix to index into, and then we go through each element in the array and fill it in with the appropriate value by calling to the coeff(int row, int col) function we defined above. Easy!

Now we can use the same setup file from the last Cython post to build our shared library:

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

setup(
  name = 'Test',
  ext_modules=[ 
    Extension("test", 
              sources=["test.pyx", "cpp_test.cpp"], # Note, you can link against a c++ library instead of including the source
              language="c++"),
    ],
  cmdclass = {'build_ext': build_ext},

)

Again reproduced here for hassle-freeness. Now we can compile and test it! At the IPython terminal:

run setup.py build_ext -i
import test
t = test.pyTest(10)
m = t.getNDArray(4,4)

And that’s it! We can see that that last command returns to us a (4x4) matrix with the first element set to -10.0101003, which is exactly what we did over in our C++ files. Now we’re all set to go ahead and generate matrices in C++, and then ship ’em out over to Python for plotting. Excellent.

Tagged , ,

N-dimensional Matrices in C++

I’m working on coding up an algorithm in C++, based on a Matlab implementation I wrote, and I’ve been using the Eigen library (http://eigen.tuxfamily.org/) for my linear algebra. It’s been relatively straightforward to use so far, but unfortunately it doesn’t (as far as I can tell) have any built in support for N-dimensional matrices. In fact, it appears that N-dimensional matrices aren’t really addressed by any of the standard C++ linear algebra libraries. So, to get my 4-dimensional matrix, I have to make an array of pointers pointing to an array of pointers pointing to my Eigen MatrixXf objects. I initialize my 4-D matrices using this function:

MatrixXf*** initialize_4dMatrix(MatrixXf*** matrix, int sizes[], float value=0.0) {
    matrix = new MatrixXf** [sizes[0]];
    for (int d1 = 0; d1 < sizes[0]; d1++) {
        matrix[d1] = new MatrixXf* [sizes[1]];
        for (int d2 = 0; d2 < sizes[1]; d2++) 
            matrix[d1][d2] = new (MatrixXf)(MatrixXf::Ones(sizes[2],sizes[3]).array() * value);
    }
    return matrix;
} 

and call it with:

int sizes[4] = {1,2,3,4};
MatrixXf ***m = initialize_4dMatrix(m, {1,2,3,4});

To then access the matrices, it’s (*m[d1][d2])(d3,d4); which isn’t bad for accessing one matrix at a time, but what if I want the jth element of the ith row of every matrix (i.e. element (i,j))? The Python access method m[:,:,i,:] (or just m[:,:,i]) would be ideal. Unfortunately, I don’t think there’s any way around writing another function to do it:

    static VectorXf* get_elements_across_matrix_array(MatrixXf** matrix, int size, int element[]) { 
        VectorXf* vector = new VectorXf(size);
        for (int i = 0; i < size; i++) 
            (*vector)[i] = (*matrix[i])(element[0],element[1]);
        return vector;
    }

And I would need to write another function that returns a MatrixXf pointer if I wanted to access the ith rows of each matrix.

To delete these matrices I also need to pass into my delete function the sizes of the array. It’s a pain that it can’t dynamically discover the array size, but I don’t think there’s any way around it in C++.

static void delete_4dMatrix(MatrixXf*** matrix, int sizes[]) {
    for (int d1 = 0; d1 < sizes[0]; d1++) {
        for (int d2 = 0; d2 < sizes[1]; d2++)
            delete matrix[d1][d2];
    }
}

Additionally, I need to have a separate function for initializing my 3D matrices, and any other dimension that I would need the way I have it set up now. Less than ideal. My guess is that there’s a elegant way to handle this with void pointers, but at the moment I only need 3D and 4D matrices so it’s not worth the effort for me to set it up.

So there are a few obstacles to using N-dimensional matrices with Eigen, but it’s really not a problem once you get a handle on it and how to properly reference different elements and matrices through your (N-2)-dimensional array of pointers.

Tagged , ,

Cython Journey Part 1: The Basics

I’ve recently begun writing up some algorithms is C++ for the speed benefits, as they require an insane amount of trials and overall simulation time, and there is a significant speedup in moving code from Matlab or Python over to C++.

BUT

It’s a large project, and I’d like to

  1. Be able to easily plot the results in Python, using matplotlib
  2. Be able to “plug and play” C++ parts with Python parts for speed comparison and error checking

So I’ve been looking into Cython, which, after some digging, seems a very viable candidate for building up this interface. I’ve been following the examples and whatnot they have on the Cython page (http://docs.cython.org/src/userguide/wrapping_CPlusPlus.html#a-simple-tutorial), but it tends not to be abundantly clear what the finished version of the files are supposed to look like. So I’ve been writing up my own as I go along based on what their examples and trial and error.

I’ve called the project test. I’ve created two files to hold the C++ code, cpp_test.h and cpp_test.cpp, which look like this:
cpp_test.h

#ifndef TEST_H
#define TEST_H

class Test { 
public: 
    int test1; 
    Test();
    Test(int test1); 
    ~Test(); 
    int returnFour(); 
    int returnFive();
    Test operator+(const Test& other); 
    Test operator-(const Test& other);
};
#endif

and cpp_test.cpp

#include "cpp_test.h"

Test::Test() { 
    test1 = 0;
}

Test::Test(int test1) { 
    this->test1 = test1; 
}

Test::~Test() { }

int Test::returnFour() { return 4; }

int Test::returnFive() { return returnFour() + 1; }

Test Test::operator+(const Test& other) { 
    return Test(test1 += other.test1);
}

Test Test::operator-(const Test& other) { 
    return Test(test1 -= other.test1);
}

To start, we need a wrapper .pyx file that will specify how we actually interface between Python and the above C++ code. So, we make a new file, and call it “test.pyx”.

  • Important! Do not name the .pyc file the same as your .cpp file! When the cythonize function is called later, a new .cpp file with the same name as your .pyx file is generated. If you name both your original .cpp and your .pyx file the same, your .cpp file will be overwritten! Which is bad.

The first thing we need to do in our test.pyx file (which is written in Cython) is redefine all of the variables and functions that we want to use from our C++ code, including the constructors. In this example we’re going to expose the operator functions and returnFive(), which calls returnFour(), but not returnFour() itself. All of our C++ interaction is defined in this .pyx file.
Here is the first part of test.pyx:

cdef extern from "cpp_test.h": 
    cdef cppclass Test:
        Test()
        Test(int test1)
        int test1
        int returnFive()
        Test add "operator+"(Test other) 
        Test sub "operator-"(Test other)

Here we have a redeclaration of everything that we want to use from our C++ code, which is everything except returnFour(). All of the class variables and functions are defined inside cdef cppclass Test. One of the first things that jumps out from looking at this code is that we can rename C++ commands with Python names here, by declaring a name, and then the mapped C++ code inside quotation marks. The above code gives Cython a handle into the C++ code, but we are not done! We need to make this functionality available to regular ol’ Python.

So! Still inside test.pyx, we need to write some more code. I’ll show the code and then walk through it.

cdef class pyTest: 
    cdef Test* thisptr # hold a C++ instance
    def __cinit__(self, int test1):
        self.thisptr = new Test(test1)
    def __dealloc__(self):
        del self.thisptr

    def __add__(pyTest left, pyTest other):
        cdef Test t = left.thisptr.add(other.thisptr[0])
        cdef pyTest tt = pyTest(t.test1)
        return tt
    def __sub__(pyTest left, pyTest other):
        cdef Test t = left.thisptr.sub(other.thisptr[0])
        cdef pyTest tt = pyTest(t.test1)
        return tt

    def __repr__(self): 
        return "pyTest[%s]" % (self.thisptr.test1)

    def returnFive(self):
        return self.thisptr.returnFive()

    def printMe(self):
        return "hello world"

Let’s look at what’s going on. We define a new class, and it’s called Test. The first thing we do is then define a pointer to the class defined above in test.pyc, c_Test. This will be how we access the variables and functions of our C++ code. The first method we define is __cinit__, which is the very first thing called when we make an instance of Test, called even before __init__. It’s important to do this here, because of some issues that come up with __init__ itself not being called soon enough for instantiating our pointer, thisptr. To create pointers to the class we simply call new Test(int test1). The next method we define is __dealloc__, which just lets us delete our pointer, which calls the deconstructor of our C++ class and frees memory appropriately, through del self.thisptr.

After these, we create a method to call returnFive(), which we appropriately also title returnFive(). All this does is call upon thisptr to access it’s own returnFive() function, which links to the C++ returnFive() function. Note that the C++ code calls upon returnFour(), which we did not define in c_Test. This just shows that not everything has to be redefined in our .pyx file, only what we want exposed for access from Python.

The next two methods are the Cython equivalent of our operator+ and operator- C++ functions. These functions are a little trickier, because the c_Test functions return pointers to new c_Test instances, because the operator functions in our C++ code in turn returns new instances of the C++ Test class. So what we have to do is create a handle to the c_Test instance returned by the self.thisptr.add and self.thisptr.sub calls, and then create a new instance of the Test class from our .pyx code, and return that.

The __repr__ function is simply an aesthetic feature, such that when the user calls, from Python, for the Test class to print we get a nice printout.

And finally, just to show that this is a Python class we’re defining, we make another method that is strictly superfluous, and unrelated to the C++ code, just to emphasize that additional functionality specific to our Python implementation can be added here.

Now we have 3 files, cpp_test.h, cpp_test.cpp, and test.pyx. We need one more, setup.py. This is file that will compile the C++ and .pyx files together and generate a usable python module. The setup file can be written several different ways, here is template I’ve been using:
setup.py

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext

setup(
  name = 'Demos',
  ext_modules=[ 
    Extension("test", 
              sources=["test.pyx", "cpp_test.cpp"], # Note, you can link against a c++ library instead of including the source
              language="c++"),
    ],
  cmdclass = {'build_ext': build_ext},

)

The important thing here is to set the sources parameter to include the relevant .pyx and .cpp files that we just wrote.
And that’s it! Now open up IPython and punch in:

run setup.py build_ext -i

and now we can call up our test module! i.e.

from test Import pyTest
k = pyTest(10)
k.returnFive()
k + k
k - k 
k.printMe()

And that’s it! We’ve successfully wrapped a basic C++ class and made it accessible from Python! This is a good first step, we’ll continue the Cython journey in later posts.

One last thing though, if you use IPython, it’s should be noted that if you change any files and recompile with the

run setup.py build_ext -i

you need to exit and then re-enter IPython before the updated code will be recognized. If you don’t know this you can get very frustrated trying to figure out why nothing is changing!

Tagged ,