2

Background

I have created a python module that wraps a c++ program using SWIG. It works just fine, but it has a pretty serious memory leak issue that I think is a result of poorly handled pointers to large map objects. I have very little experience with c++, and I have questions as to whether delete[] can be used on an object created with new in a different function or method.

The program was written in 2007, so excuse the lack of useful c++11 tricks.


The swig extension basically just wraps a single c++ class (Matrix) and a few functions.

Matrix.h

#ifndef __MATRIX__
#define __MATRIX__

#include <string>
#include <vector>
#include <map>
#include <cmath>
#include <fstream>
#include <cstdlib>
#include <stdio.h>
#include <unistd.h>

#include "FileException.h"
#include "ParseException.h"

#define ROUND_TO_INT(n) ((long long)floor(n))
#define MIN(a,b) ((a)<(b)?(a):(b))
#define MAX(a,b) ((a)>(b)?(a):(b))

using namespace std;

class Matrix {

private:


  /**
  * Split a string following delimiters
   */
  void tokenize(const string& str, vector<string>& tokens, const string& delimiters) {

    // Skip delimiters at beginning.
    string::size_type lastPos = str.find_first_not_of(delimiters, 0);
    // Find first "non-delimiter".
    string::size_type pos     = str.find_first_of(delimiters, lastPos);

    while (string::npos != pos || string::npos != lastPos)
    {
      // Found a token, add it to the vector.
      tokens.push_back(str.substr(lastPos, pos - lastPos));
      // Skip delimiters.  Note the "not_of"
      lastPos = str.find_first_not_of(delimiters, pos);
      // Find next "non-delimiter"
      pos = str.find_first_of(delimiters, lastPos);
    }
  }


public:


  // used for efficiency tests
  long long totalMapSize;
  long long totalOp;

  double ** mat; // the matrix as it is stored in the matrix file
  int length;
  double granularity; // the real granularity used, greater than 1
  long long ** matInt; // the discrete matrix with offset
  double errorMax;
  long long *offsets; // offset of each column
  long long offset; // sum of offsets
  long long *minScoreColumn; // min discrete score at each column
  long long *maxScoreColumn; // max discrete score at each column
  long long *sum;
  long long minScore;  // min total discrete score (normally 0)
  long long maxScore;  // max total discrete score
  long long scoreRange;  // score range = max - min + 1
  long long *bestScore;
  long long *worstScore;
  double background[4];

  Matrix() {
    granularity = 1.0;
    offset = 0;
    background[0] = background[1] = background[2] = background[3] = 0.25;
  }

  Matrix(double pA, double pC, double pG, double pT) {
    granularity = 1.0;
    offset = 0;
    background[0] = pA;
    background[1] = pC;
    background[2] = pG;
    background[3] = pT;  
  }

  ~Matrix() {
      for (int k = 0; k < 4; k++ ) {
        delete[] matInt[k];
      }
      delete[] matInt;
      delete[] mat;
      delete[] offsets;
      delete[] minScoreColumn;
      delete[] maxScoreColumn;
      delete[] sum;
      delete[] bestScore;
      delete[] worstScore;
  }


  void toLogOddRatio () {
    for (int p = 0; p < length; p++) {
      double sum = mat[0][p] + mat[1][p] + mat[2][p] + mat[3][p];
      for (int k = 0; k < 4; k++) {
        mat[k][p] = log((mat[k][p] + 0.25) /(sum + 1)) - log (background[k]); 
      }
    }
  }

  void toLog2OddRatio () {
    for (int p = 0; p < length; p++) {
      double sum = mat[0][p] + mat[1][p] + mat[2][p] + mat[3][p];
      for (int k = 0; k < 4; k++) {
        mat[k][p] = log2((mat[k][p] + 0.25) /(sum + 1)) - log2 (background[k]); 
      }
    }
  }

  /**
    * Transforms the initial matrix into an integer and offseted matrix.
   */
  void computesIntegerMatrix (double granularity, bool sortColumns = true);

  // computes the complete score distribution between score min and max
  void showDistrib (long long min, long long max) {
    map<long long, double> *nbocc = calcDistribWithMapMinMax(min,max); 
    map<long long, double>::iterator iter;

    // computes p values and stores them in nbocc[length] 
    double sum = 0;
    map<long long, double>::reverse_iterator riter = nbocc[length-1].rbegin();
    while (riter != nbocc[length-1].rend()) {
      sum += riter->second;
      nbocc[length][riter->first] = sum;
      riter++;      
    }

    iter = nbocc[length].begin();
    while (iter != nbocc[length].end() && iter->first <= max) {
      //cout << (((iter->first)-offset)/granularity) << " " << (iter->second) << " " << nbocc[length-1][iter->first] << endl;
      iter ++;
    }
  }

  /**
    * Computes the pvalue associated with the threshold score requestedScore.
    */
  void lookForPvalue (long long requestedScore, long long min, long long max, double *pmin, double *pmax);

  /**
    * Computes the score associated with the pvalue requestedPvalue.
    */
  long long lookForScore (long long min, long long max, double requestedPvalue, double *rpv, double *rppv);

  /** 
    * Computes the distribution of scores between score min and max as the DP algrithm proceeds 
    * but instead of using a table we use a map to avoid computations for scores that cannot be reached
    */
  map<long long, double> *calcDistribWithMapMinMax (long long min, long long max); 

  void readMatrix (string matrix) {

    vector<string> str;
    tokenize(matrix, str, " \t|");
    this->length = 0;
    this->length = str.size() / 4;
    mat = new double*[4];
    int idx = 0;

    for (int j = 0; j < 4; j++) {
      this->mat[j] = new double[this->length];
      for (int i = 0; i < this->length; i++) {
        mat[j][i] = atof(str.at(idx).data());
        idx++;
      }
    }

    str.clear();

  }

}; /* Matrix */

#endif

Matrix.cpp

#include "Matrix.h"

#define MEMORYCOUNT

void Matrix::computesIntegerMatrix (double granularity, bool sortColumns) {
  double minS = 0, maxS = 0;
  double scoreRange;

  // computes precision
  for (int i = 0; i < length; i++) {
    double min = mat[0][i];
    double max = min;
    for (int k = 1; k < 4; k++ )  {
      min = ((min < mat[k][i])?min:(mat[k][i]));
      max = ((max > mat[k][i])?max:(mat[k][i]));
    }
    minS += min;
    maxS += max;
  } 

  // score range
  scoreRange = maxS - minS + 1;

  if (granularity > 1.0) {
    this->granularity = granularity / scoreRange;
  } else if (granularity < 1.0) {
    this->granularity = 1.0 / granularity;
  } else {
    this->granularity = 1.0;
  }

  matInt = new long long *[length];
  for (int k = 0; k < 4; k++ ) {
    matInt[k] = new long long[length];
    for (int p = 0 ; p < length; p++) {
      matInt[k][p] = ROUND_TO_INT((double)(mat[k][p]*this->granularity)); 
    }
  }

  this->errorMax = 0.0;
  for (int i = 1; i < length; i++) {
    double maxE = mat[0][i] * this->granularity - (matInt[0][i]);
    for (int k = 1; k < 4; k++) {
      maxE = ((maxE < mat[k][i] * this->granularity - matInt[k][i])?(mat[k][i] * this->granularity - (matInt[k][i])):(maxE));
    }
    this->errorMax += maxE;
  }

  if (sortColumns) {
    // sort the columns : the first column is the one with the greatest value
    long long min = 0;
    for (int i = 0; i < length; i++) {
      for (int k = 0; k < 4; k++) {
        min = MIN(min,matInt[k][i]);
      }
    }
    min --;
    long long *maxs = new long long [length];
    for (int i = 0; i < length; i++) {
      maxs[i] = matInt[0][i];
      for (int k = 1; k < 4; k++) {
        if (maxs[i] < matInt[k][i]) {
          maxs[i] = matInt[k][i];
        }
      }
    }
    long long **mattemp = new long long *[4];
    for (int k = 0; k < 4; k++) {        
      mattemp[k] = new long long [length];
    }
    for (int i = 0; i < length; i++) {
      long long max = maxs[0];
      int p = 0;
      for (int j = 1; j < length; j++) {
        if (max < maxs[j]) {
          max = maxs[j];
          p = j;
        }
      }
      maxs[p] = min;
      for (int k = 0; k < 4; k++) {        
        mattemp[k][i] = matInt[k][p];
      }
    }

    for (int k = 0; k < 4; k++)  {
      for (int i = 0; i < length; i++) {
        matInt[k][i] = mattemp[k][i];
      }
    }

    for (int k = 0; k < 4; k++) {        
      delete[] mattemp[k];
    }
    delete[] mattemp;
    delete[] maxs;
  }

  // computes offsets
  this->offset = 0;
  offsets = new long long [length];
  for (int i = 0; i < length; i++) {
    long long min = matInt[0][i];
    for (int k = 1; k < 4; k++ )  {
      min = ((min < matInt[k][i])?min:(matInt[k][i]));
    }
    offsets[i] = -min;
    for (int k = 0; k < 4; k++ )  {
      matInt[k][i] += offsets[i];  
    }
    this->offset += offsets[i];
  }

  // look for the minimum score of the matrix for each column
  minScoreColumn = new long long [length];
  maxScoreColumn = new long long [length];
  sum            = new long long [length];
  minScore = 0;
  maxScore = 0;
  for (int i = 0; i < length; i++) {
    minScoreColumn[i] = matInt[0][i];
    maxScoreColumn[i] = matInt[0][i];
    sum[i] = 0;
    for (int k = 1; k < 4; k++ )  {
      sum[i] = sum[i] + matInt[k][i];
      if (minScoreColumn[i] > matInt[k][i]) {
        minScoreColumn[i] = matInt[k][i];
      }
      if (maxScoreColumn[i] < matInt[k][i]) {
        maxScoreColumn[i] = matInt[k][i];
      }
    }
    minScore = minScore + minScoreColumn[i];
    maxScore = maxScore + maxScoreColumn[i];
    //cout << "minScoreColumn[" << i << "] = " << minScoreColumn[i] << endl;
    //cout << "maxScoreColumn[" << i << "] = " << maxScoreColumn[i] << endl;
  }
  this->scoreRange = maxScore - minScore + 1;

  bestScore = new long long[length];
  worstScore = new long long[length];
  bestScore[length-1] = maxScore;
  worstScore[length-1] = minScore;
  for (int i = length - 2; i >= 0; i--) {
    bestScore[i]  = bestScore[i+1]  - maxScoreColumn[i+1];
    worstScore[i] = worstScore[i+1] - minScoreColumn[i+1];
  }


}




/**
* Computes the pvalue associated with the threshold score requestedScore.
 */
void Matrix::lookForPvalue (long long requestedScore, long long min, long long max, double *pmin, double *pmax) {

  map<long long, double> *nbocc = calcDistribWithMapMinMax(min,max); 
  map<long long, double>::iterator iter;


  // computes p values and stores them in nbocc[length] 
  double sum = nbocc[length][max+1];
  long long s = max + 1;
  map<long long, double>::reverse_iterator riter = nbocc[length-1].rbegin();
  while (riter != nbocc[length-1].rend()) {
    sum += riter->second;
    if (riter->first >= requestedScore) s = riter->first;
    nbocc[length][riter->first] = sum;
    riter++;      
  }
  //cout << "   s found : " << s << endl;

  iter = nbocc[length].find(s);
  while (iter != nbocc[length].begin() && iter->first >= s - errorMax) {
    iter--;      
  }
  //cout << "   s - E found : " << iter->first << endl;

#ifdef MEMORYCOUNT
  // for tests, store the number of memory bloc necessary
  for (int pos = 0; pos <= length; pos++) {
    totalMapSize += nbocc[pos].size();
  }
#endif

  *pmax = nbocc[length][s];
  *pmin = iter->second;

}



/**
* Computes the score associated with the pvalue requestedPvalue.
 */
long long Matrix::lookForScore (long long min, long long max, double requestedPvalue, double *rpv, double *rppv) {

  map<long long, double> *nbocc = calcDistribWithMapMinMax(min,max); 
  map<long long, double>::iterator iter;

  // computes p values and stores them in nbocc[length] 
  double sum = 0.0;
  map<long long, double>::reverse_iterator riter = nbocc[length-1].rbegin();
  long long alpha = riter->first+1;
  long long alpha_E = alpha;
  nbocc[length][alpha] = 0.0;
  while (riter != nbocc[length-1].rend()) {
    sum += riter->second;
    nbocc[length][riter->first] = sum;
    if (sum >= requestedPvalue) { 
      break;
    }
    riter++;      
  }
  if (sum > requestedPvalue) {
    alpha_E = riter->first;
    riter--;
    alpha = riter->first; 
  } else {
    if (riter == nbocc[length-1].rend()) { // path following the remark of the mail
      riter--;
      alpha = alpha_E = riter->first;
    } else {
      alpha = riter->first;
      riter++;
      sum += riter->second;
      alpha_E = riter->first;
    }
    nbocc[length][alpha_E] = sum;  
    //cout << "Pv(S) " << riter->first << " " << sum << endl;   
  } 

#ifdef MEMORYCOUNT
  // for tests, store the number of memory bloc necessary
  for (int pos = 0; pos <= length; pos++) {
    totalMapSize += nbocc[pos].size();
  }
#endif

  if (alpha - alpha_E > errorMax) alpha_E = alpha;

  *rpv = nbocc[length][alpha];
  *rppv = nbocc[length][alpha_E];   

  delete[] nbocc;
  return alpha;

}


// computes the distribution of scores between score min and max as the DP algrithm proceeds 
// but instead of using a table we use a map to avoid computations for scores that cannot be reached
map<long long, double> *Matrix::calcDistribWithMapMinMax (long long min, long long max) { 

  // maps for each step of the computation
  // nbocc[length] stores the pvalue
  // nbocc[pos] for pos < length stores the qvalue
  map<long long, double> *nbocc = new map<long long, double> [length+1];
  map<long long, double>::iterator iter;

  long long *maxs = new long long[length+1]; // @ pos i maximum score reachable with the suffix matrix from i to length-1

  maxs[length] = 0;
  for (int i = length-1; i >= 0; i--) {
    maxs[i] = maxs[i+1] + maxScoreColumn[i];
  }

  // initializes the map at position 0
  for (int k = 0; k < 4; k++) {
    if (matInt[k][0]+maxs[1] >= min) {
      nbocc[0][matInt[k][0]] += background[k];
    }
  }

  // computes q values for scores greater or equal than min
  nbocc[length-1][max+1] = 0.0;
  for (int pos = 1; pos < length; pos++) {
    iter = nbocc[pos-1].begin();
    while (iter != nbocc[pos-1].end()) {
      for (int k = 0; k < 4; k++) {
        long long sc = iter->first + matInt[k][pos];
        if (sc+maxs[pos+1] >= min) {
          // the score min can be reached
          if (sc > max) {
            // the score will be greater than max for all suffixes
            nbocc[length-1][max+1] += nbocc[pos-1][iter->first] * background[k]; //pow(4,length-pos-1) ;
            totalOp++;
          } else {              
            nbocc[pos][sc] += nbocc[pos-1][iter->first] * background[k];
            totalOp++;
          }
        } 
      }
      iter++;      
    }      
    //cerr << "        map size for " << pos << " " << nbocc[pos].size() << endl;
  }

  delete[] maxs;

  return nbocc;


}

pytfmpval.i

%module pytfmpval
%{
#include "../src/Matrix.h"
#define SWIG_FILE_WITH_INIT
%}

%include "cpointer.i"
%include "std_string.i"
%include "std_vector.i"
%include "typemaps.i"
%include "../src/Matrix.h"

%pointer_class(double, doublep)
%pointer_class(int, intp)

%nodefaultdtor Matrix;

The c++ functions are called in a python module.


I worry that nbocc in Matrix.cpp is not being properly dereferenced or deleted. Is this use valid?

I have tried using gc.collect() and I am using the multiprocessing module as recommended in this question to call these functions from my python program. I've also tried deleting the Matrix object from within python to no avail.

I'm out of characters, but will provide any additional needed info in the comments as well as I can.

UPDATE: I've removed all of the python code, as it wasn't the issue and made for an absurdly long post. As I stated in the comments below, this was ultimately solved by taking the suggestion of many users and creating a minimal example that exhibited the issue in pure C++. I then used valgrind to identify the problematic pointers created with new and made sure that they were properly dereferenced. This fixed almost all memory leaks. One remains, but it leaks only a few hundred bytes over thousands of iterations and would require refactoring the entire Matrix class, which simply isn't worth the time for what it is. Bad practice, I know. To any other newbie in C++ out there, seriously try to avoid dynamic memory allocation or utilize std::unique_ptr or std::shared_ptr.

Thanks again to everyone who provided input and suggestions.

Jared Andrews
  • 272
  • 2
  • 12
  • 5
    You need to show a *minimal* complete example we can try with a debugger on our own computers that reproduces the issue and not random incomplete fragments of a large fragment. – Flexo Jul 07 '17 at 16:43
  • @Flexo Sorry, the complete setup is pretty extensive and likely overkill. I've updated my question with the source files for the `Matrix` class and how they're being utilized by the python module. I've also spent a lot of time reading about pointers and dynamic memory allocation and updated my question with related concerns. Thanks for your time. – Jared Andrews Jul 07 '17 at 19:13
  • I am not sure how your `Matrix.cpp` compiles without error. In `Matrix::lookForPvalue` you do `map nbocc = calcDistribWithMapMinMax(min,max);` which should cause compilation error since LHS `nbocc` is not a type of `map *` which is what `calcDistribWithMapMinMax` returns. – Griffin Jul 11 '17 at 15:24
  • That's a typo, my mistake, I will edit the code. – Jared Andrews Jul 11 '17 at 15:26
  • @JaredAndrews I am no expert in SWIG, but, could you try removing `%nodefaultdtor Matrix;` and see if you can still reproduce the memory leak? – Griffin Jul 11 '17 at 15:29
  • @Griffin Yes, it will still be there. I added that after trying to explicitely define a destructor, but the problem existed prior to that. None of the `Matrix` attributes themselves consume much memory and appear properly cleaned up on the Python side. – Jared Andrews Jul 11 '17 at 15:33
  • 1
    @JaredAndrews How have you checked memory is properly cleaned up on Pythod side? Can you try tracing the leak with [this suggested tool](http://tech.labs.oliverwyman.com/blog/2008/11/14/tracing-python-memory-leaks/)? – Griffin Jul 11 '17 at 15:43
  • @Griffin. This looks useful and I'll test it out. I used [Pympler](https://pythonhosted.org/Pympler/) to check out the size of objects while debugging. I will try tracing it will this tool as well. – Jared Andrews Jul 11 '17 at 15:48
  • @Griffin So using both that tool and Pympler yield basically the same result. `Matrix` objects themselves have a very small memory footprint and are deleted properly, *at least* as can be seen from the python side. [This Pympler image](http://i.imgur.com/jLyXwwU.png) shows all `Matrix` objects right after instantiation and then deletion. I think it's going to be difficult to assay the leak from the python side, so I'm trying to get Valgrind to play nice and see what it yields. – Jared Andrews Jul 11 '17 at 17:22
  • @JaredAndrews could you share what you have with Valgrind found with us? – Griffin Jul 12 '17 at 09:18
  • 1
    This can/ **should** be debugged outside _Python_. Simply add a _main.cpp_ file to the "project", that instantiates the `Matrix` class and call its methods. The rule is simple: for every `new` call there should be one `delete[]` counterpart (brackets might not be necessary). E.g. all matrix type (double pointer - `**`) variables that are allocated must be freed (like `matInt` from destructor). `mattemp` on the other hand, produces such a memory leak (there might be more). – CristiFati Jul 13 '17 at 00:03
  • A general comment: you could avoid the use of `new`/`delete` here, and use `std::vector` or `std::array` instead, then you are likely to avoid memory leaks. – Mine Jul 13 '17 at 07:17
  • @Griffin Valgrind confirmed the memory leak, but also yielded a ton of noise even with the python suppressions. I'm creating a purely c++ example as suggested so that I can test the solution in Some Who Call Me Tim's answer below. – Jared Andrews Jul 13 '17 at 16:19
  • You say it was written in 2007 so doesn’t have C++11 tricks, but is there a reason you couldn’t write it to use those now? You could get rid of almost all the explicit `new`s and `delete`s if you did that, and the fewer of those you have the fewer chances there are for memory leaks to sneak in. – Daniel H Jul 13 '17 at 22:44

3 Answers3

2

Two questions are in play here: managing memory in C++, and then nudging the C++ side from the Python side to clean up. I'm guessing SWIG is generating a wrapper for the Matrix destructor and calling the destructor at some useful time. (I might convince myself of that by having the dtor make some noise.) That should handle the second question.

So let's focus on the C++ side. Passing around a bare map * is a well-known invitation to mischief. Here are two alternatives.

Alternative one: make the map a member of Matrix. Then it gets cleaned up automatically by ~Matrix(). This is the easiest thing. If the lifetime of the map does not exceed the lifetime of the Matrix, then this route will work.

Alternative two: if the map needs to persist after the Matrix object, then instead of passing around map *, use a shared pointer, std::shared_ptr<map>. The shared pointer reference counts the pointee (i.e. the dynamically allocated Matrix). When the ref count goes to zero, it deletes the underlying object.

They both build on the rule to allocate resources (memory in this case) in constructors and deallocate in destructors. This is called RAII (Resource Allocation Is Initialization). Another application of RAII in your code would be to use std::vector<long long> offsets instead of long long *offsets etc. Then you just resize the vectors as needed. When the Matrix is destroyed, the vectors are deleted with no intervention on your part. For the matrix, you could use a vector of vectors, and so on.

  • Thanks. Yes, SWIG does just that, and at this point, I'm quite certain the `Matrix` class is being cleaned up properly. I had contemplated your first approach and think it's likely the best solution. I had come across shared pointers previously, but this isn't c++11, so it's not an option unless I use something like `boost`, which I'm less inclined to do. The vector usage is also a good idea and something I will try to implement. – Jared Andrews Jul 13 '17 at 16:13
2

It’s hard to follow what’s happening, but I’m pretty sure your matrices are not being cleaned up correctly.

In readMatrix, you have a loop over j which contains the line this->mat[j] = new double[this->length];. This allocates memory, which mat[j] points to. This memory needs to be freed at some point, by calling delete[] mat[j] (or some other loop variable). However, in the destructor, you just call delete[] mat, which leaks all of the arrays inside it.

Some general suggestions on cleaning this up:

  1. If you know the bounds of an array, such as that matInt will always have a length of 4, you should declare it with that fixed length (long long* matInt[4] will make an array of four pointers to long long, each of which could be a pointer to an array); this will mean you don’t need to either new or delete it.
  2. If you have a double pointer like double ** mat, and you allocate both the first and second layers of pointers with new[], you need to deallocate the inner layer with delete[] (and you need to do it before you delete[] the outer layer).
  3. If you still have trouble, more of your code will be clear if you remove the methods which don’t seem relevant to the problem. For example, toLogOddRatio doesn’t allocate or deallocate memory at all; it almost certainly isn’t contributing to the problem and you can remove it from the code you post here (once you’ve removed the parts which you think don’t contribute, test again to make sure the problem’s still there; if not then you know that it was one of those parts somehow causing the leak).
Daniel H
  • 7,223
  • 2
  • 26
  • 41
  • 1
    Though all of these comments and answers were helpful, this one was the one that ultimately best addressed the issues I was facing. As previously suggested by @Flexo and others, I created a minimal example in c++ that still had the issue and used `valgrind` to find where the leaks were. This allowed me to systematically remove/fix all `new` and `delete` statements that we causing the leaks. Thanks again to everyone for your help. The source code for the corrected program can be seen here: https://github.com/j-andrews7/pytfmpval – Jared Andrews Sep 01 '17 at 19:12
  • @JaredAndrews I’m glad I could help! One tip which I forgot to mention here is that it's often easier and less error-prone to use C++ classes like [`std::unique_ptr`](http://en.cppreference.com/w/cpp/memory/unique_ptr) instead of trying to do memory management yourself; unfortunately, I don’t know if those play nicely with SWIG. – Daniel H Sep 01 '17 at 19:30
  • 1
    Yeah, they're supported (supposedly, though I read plenty of threads that had issues) but I really wasn't looking forward to reworking everything with the different pointers - this was my first contact with C++ and this was meant to be a quick task that turned into something that was much more of a headache. If I ever wrote anything from scratch, I'd be sure to utilize the features of C++11 though. Thanks again! – Jared Andrews Sep 01 '17 at 19:34
1

to answer your question, yes you can use delete on diffrent function or method. and you should, any memory you allocate in c/c++ you need to free (delete in c++ lingo)

python isn't aware of this memory, it's not a python object, so gc.collect() won't help. you should add a c function that would take a Matrix struct and free/delete the memory use on that struct. and call it from python, swig in not handling memory allocation (only for the objects swig creates)

I would recommended looking into newer packages other then swig, like cython or cffi (or even NumPy matrix handling, I've heard he's good at)

Fruch
  • 408
  • 5
  • 18
  • Thanks for the first answer. I actually did just what you suggested, but the `Matrix` object isn't actually the issue, but rather the extensive use of pointers to maps of maps and issues dereferencing them properly. I looked at other options (Boost, Cython, ctypes among them), but Swig seemed like the most straightforward one. – Jared Andrews Jul 11 '17 at 15:24