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.