The main problem here was that you were getting copies of nump
arrays(numpy_values, numpy_sd, numpy_relax_times, numpy_scaling_matrix)
that you own from PyArray_ContiguousFromObject and never decrefing
them, so their reference counts never reached zero... However, if you
Py_DECREF these objects at the end of setup they will be collected and
the double arrays values, sd, relax_time, scaling_matrix which you
still have references to will no longer be allocated to these objects
and will contain junk if the memory allocator has reallocated stuff
where they were. So I have set it up so that you keep both
numpy_values, numpy_sd, numpy_relax_times, numpy_scaling_matrix and
values, sd, relax_time, scaling_matrix in global varaibles and only
call Py_DECREF at the start of setup to junk them just before you
reallocate them (actually I call Py_XDECREF as they will be null the
first time the setup function is called).
So here are the patched files c math_fns which cure the memeory leaks Ed.
The changes I have made are as follows
1. for c_chi2 and exponential I have made them both into proper
functions with no references to any globals
2. for the relax_fit.c file I have basically done the following
a. I have removed the global visibilty of numpy_params and made sure
that i call Py_DECREF in each function where you call
PyArray_ContiguousFromObject
b. I have made sure that there are numpy_values, numpy_sd,
numpy_relax_times, numpy_scaling_matrix.
c. I call Py_XDECREF for each of numpy_values, numpy_sd,
numpy_relax_times, numpy_scaling_matrix so that we don't orphan them
when we call PyArray_ContiguousFromObject but ignore null objects the
first time we are called.
3. note dfunc has a rather funny return ;-)
/* Test code (convert aaa to a Numeric array */
/* aaa_numpy = (PyArrayObject *) PyArray_FromDimsAndData(1,
num_params, PyArray_DOUBLE, aaa_pointer); */
/*aaa_numpy = (PyArrayObject *) PyArray_FromDims(1, &num_params,
PyArray_DOUBLE);
aaa_pointer = (double *) aaa_numpy->data;*/
/* Fill the Numeric array */
/*for (i = 0; i < 2; i++)
aaa_pointer[i] = aaa[i];*/
Py_DECREF(numpy_params);
return NULL;
return PyArray_Return(aaa_numpy);
which may cause some more careful compilers to choke (but not g77 ;-)
4. Other comments (more things i could also sort at some point)
a this should really be an object, as it has state and a lifetime.
If I get time I will convert it...
b. MAXTIMES shouldn't really exist, data arrays should be allocated
on the fly I guess
regards
gary
--
-------------------------------------------------------------------
Dr Gary Thompson
Astbury Centre for Structural Molecular Biology,
University of Leeds, Astbury Building,
Leeds, LS2 9JT, West-Yorkshire, UK Tel. +44-113-3433024
email: garyt@xxxxxxxxxxxxxxx Fax +44-113-2331407
-------------------------------------------------------------------
/*
* Copyright (C) 2003, 2006 Edward d'Auvergne
*
* This file is part of the program relax.
*
* relax is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* relax is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with relax; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
#include <stdio.h>
#include <math.h>
#define square(x) (x)*(x)
double chi2(double *values, double *sd, double *back_calc, int num_times) {
/* Function to calculate the chi-squared value.
The chi-sqared equation
~~~~~~~~~~~~~~~~~~~~~~~
_n_
\ (yi - yi()) ** 2
Chi2() = > ----------------
/__ sigma_i ** 2
i=1
where:
yi are the values of the measured data set.
yi() are the values of the back calculated data set.
sigma_i are the values of the error set.
The chi-squared value is returned.
*/
int i;
double chi2 = 0.0;
/* Loop over the time points and sum the chi-squared components */
for (i = 0; i < num_times; ++i) {
chi2 = chi2 + square((values[i] - back_calc[i]) / sd[i]);
}
return chi2;
}
/*
* Copyright (C) 2006 Edward d'Auvergne
*
* This file is part of the program relax.
*
* relax is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* relax is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with relax; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
#ifndef RELAX_C_CHI2
#define RELAX_C_CHI2
double chi2(double *values, double *sd, double *back_calc, int num_times);
#endif
/*
* Copyright (C) 2006 Edward d'Auvergne
*
* This file is part of the program relax.
*
* relax is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* relax is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with relax; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
/* This include must come first */
#include <Python.h>
/* The Numeric array object header file, must come second */
#include <Numeric/arrayobject.h>
/* The header for all functions which will be called */
#include "relax_fit.h"
void exponential(double *params, double *relax_times, double *back_calc, int
num_times) {
/* Function to back calculate the peak intensities.
*/
/* Declarations */
double Rx, I0;
int i;
/* Loop over the time points */
/* for (i = 0; i < num_times; i++) { */
for (i = 0; i < num_times; i++) {
/* Zero Rx value */
if (params[0] == 0.0)
back_calc[i] = 0.0;
/* Back calculate */
else
back_calc[i] = params[1] * exp(-relax_times[i] * params[0]);
}
}
/*
* Copyright (C) 2006 Edward d'Auvergne
*
* This file is part of the program relax.
*
* relax is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* relax is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with relax; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
#ifndef RELAX_EXPONENTIAL
#define RELAX_EXPONENTIAL
void exponential(double *params, double *relax_times, double *back_calc, int
num_times);
#endif
/*
* Copyright (C) 2006 Edward d'Auvergne
*
* This file is part of the program relax.
*
* relax is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* relax is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with relax; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
/* This include must come first */
#include <Python.h>
/* The Numeric array object header file, must come second */
#include <Numeric/arrayobject.h>
/* The header for all functions which will be called */
#include "relax_fit.h"
/* functions for chi2 and exponential */
#include "c_chi2.h"
#include "exponential.h"
static PyObject *
setup(PyObject *self, PyObject *args, PyObject *keywords) {
/* Python declarations */
PyObject *values_arg, *sd_arg, *relax_times_arg, *scaling_matrix_arg;
extern PyArrayObject *numpy_values, *numpy_sd, *numpy_relax_times,
*numpy_scaling_matrix;
/* Normal declarations */
extern double *values, *sd, *relax_times, *scaling_matrix;
extern double relax_time_array;
extern int num_params, num_times;
/* The keyword list */
static char *keyword_list[] = {"num_params", "num_times", "values", "sd",
"relax_times", "scaling_matrix", NULL};
/* Parse the function arguments */
if (!PyArg_ParseTupleAndKeywords(args, keywords, "iiOOOO", keyword_list,
&num_params, &num_times, &values_arg, &sd_arg, &relax_times_arg,
&scaling_matrix_arg))
return NULL;
Py_XDECREF(numpy_values);
Py_XDECREF(numpy_sd);
Py_XDECREF(numpy_relax_times);
Py_XDECREF(numpy_scaling_matrix);
/* Make the Numeric arrays contiguous */
numpy_values = (PyArrayObject *) PyArray_ContiguousFromObject(values_arg,
PyArray_DOUBLE, 1, 1);
numpy_sd = (PyArrayObject *) PyArray_ContiguousFromObject(sd_arg,
PyArray_DOUBLE, 1, 1);
numpy_relax_times = (PyArrayObject *)
PyArray_ContiguousFromObject(relax_times_arg, PyArray_DOUBLE, 1, 1);
numpy_scaling_matrix = (PyArrayObject *)
PyArray_ContiguousFromObject(scaling_matrix_arg, PyArray_DOUBLE, 2, 2);
/* Pointers to the Numeric arrays */
values = (double *) numpy_values->data;
sd = (double *) numpy_sd->data;
relax_times = (double *) numpy_relax_times->data;
scaling_matrix = (double *) numpy_scaling_matrix->data;
/* Return nothing */
Py_INCREF(Py_None);
return Py_None;
}
static PyObject *
func(PyObject *self, PyObject *args) {
/* Function for calculating and returning the chi-squared value.
*
* Firstly the back calculated intensities are generated, then the
chi-squared statistic is
* calculated
*/
/* Declarations */
PyObject *arg1;
PyArrayObject *numpy_params;
double* params;
/* Parse the function arguments, the only argument should be the
parameter array */
if (!PyArg_ParseTuple(args, "O", &arg1))
return NULL;
/* Convert the Numeric array to be contiguous */
numpy_params = (PyArrayObject *) PyArray_ContiguousFromObject(arg1,
PyArray_DOUBLE, 1, 1);
/* Pointers to the Numeric arrays */
params = (double *) numpy_params->data;
/* Back calculated the peak intensities */
exponential(params, relax_times, back_calc, num_times);
Py_DECREF(numpy_params);
/* Calculate and return the chi-squared value */
return Py_BuildValue("f", chi2(values,sd,back_calc,num_times));
}
static PyObject *
dfunc(PyObject *self, PyObject *args) {
/* Function for calculating and returning the chi-squared gradient. */
/* Declarations */
PyObject *arg1;
PyArrayObject *numpy_params;
/* Temp Declarations */
PyArrayObject *aaa_numpy;
double aaa[MAXPARAMS] = {1.0, 2.0};
double *aaa_pointer;
int i;
double* params;
/* Parse the function arguments, the only argument should be the
parameter array */
if (!PyArg_ParseTuple(args, "O", &arg1))
return NULL;
/* Convert the Numeric array to be contiguous */
numpy_params = (PyArrayObject *) PyArray_ContiguousFromObject(arg1,
PyArray_DOUBLE, 1, 1);
/* Pointers to the Numeric arrays */
params = (double *) numpy_params->data;
/* Back calculated the peak intensities */
exponential(params, relax_times, back_calc, num_times);
/* Test code (convert aaa to a Numeric array */
/* aaa_numpy = (PyArrayObject *) PyArray_FromDimsAndData(1, num_params,
PyArray_DOUBLE, aaa_pointer); */
/*aaa_numpy = (PyArrayObject *) PyArray_FromDims(1, &num_params,
PyArray_DOUBLE);
aaa_pointer = (double *) aaa_numpy->data;*/
/* Fill the Numeric array */
/*for (i = 0; i < 2; i++)
aaa_pointer[i] = aaa[i];*/
Py_DECREF(numpy_params);
return NULL;
return PyArray_Return(aaa_numpy);
}
static PyObject *
d2func(PyObject *self, PyObject *args) {
/* Function for calculating and returning the chi-squared Hessian. */
return Py_BuildValue("f", 0.0);
}
static PyObject *
back_calc_I(PyObject *self, PyObject *args) {
/* Function for returning as a Numeric array the back calculated peak
intensities */
/* Declarations */
extern double back_calc[];
extern int num_times;
int i;
PyObject *back_calc_py = PyList_New(num_times);
assert(PyList_Check(back_calc_py));
/* Copy the values out of the C array into the Python array */
for (i = 0; i < num_times; i++)
PyList_SetItem(back_calc_py, i, Py_BuildValue("f", back_calc[i]));
/* Return the Numeric array */
return back_calc_py;
}
/* The method table for the functions called by Python */
static PyMethodDef relax_fit_methods[] = {
{"setup", (PyCFunction)setup, METH_VARARGS | METH_KEYWORDS, "The main
relaxation curve fitting setup function."},
{"func", func, METH_VARARGS},
{"dfunc", dfunc, METH_VARARGS},
{"d2func", d2func, METH_VARARGS},
{"back_calc_I", back_calc_I, METH_VARARGS},
{NULL, NULL, 0, NULL} /* Sentinel */
};
/* Initialise as a Python module */
PyMODINIT_FUNC
initrelax_fit(void)
{
(void) Py_InitModule("relax_fit", relax_fit_methods);
/* Import the Numeric array module. This is essential. */
import_array();
}
/*
* Copyright (C) 2006 Edward d'Auvergne
*
* This file is part of the program relax.
*
* relax is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* relax is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with relax; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
/* Required for the Python/C API??? */
#define PY_ARRAY_UNIQUE_SYMBOL numarray
/* The maximum number of parameters for this function */
#define MAXPARAMS 3
/* The maximum number of spectral time points */
#define MAXTIMES 30
/****************************************/
/* External, hence permanent, variables */
/****************************************/
/* Variables sent to the setup function to be stored for later use */
PyArrayObject *numpy_values, *numpy_sd, *numpy_relax_times,
*numpy_scaling_matrix;
int num_params, num_times;
double *sd;
/* Variables sent to 'func', 'dfunc', and 'd2func' during optimisation */
/*PyArrayObject *numpy_params;*/
/* Pointers to contiguous PyArrayObjects */
double *values, *sd, *relax_times, *scaling_matrix;
/*double *params;*/
/* Variables used for storage during the function calls of optimisation */
double back_calc[MAXTIMES];
math_fns c code