Author: bugman
Date: Mon Aug 25 15:28:43 2014
New Revision: 25249
URL: http://svn.gna.org/viewcvs/relax?rev=25249&view=rev
Log:
Implementation of the target_functions.relax_fit.jacobian() function.
This follows from the discussions at
http://thread.gmane.org/gmane.science.nmr.relax.devel/6807.
The function will calculate the Jacobian matrix for the exponential
curve-fitting module. The
Jacobian can be used to directly calculate the covariance matrix, for example
as described at
https://www.gnu.org/software/gsl/manual/html_node/Computing-the-covariance-matrix-of-best-fit-parameters.html.
The Jacobian is calculated using the help of the new exponential_dI() and
exponential_dR() functions
in the target_functions/exponential.c file. These calculate the partial
derivatives of the
exponential curve with respect to each model parameter separately.
The implementation still needs testing and debugging.
Modified:
trunk/target_functions/exponential.c
trunk/target_functions/exponential.h
trunk/target_functions/relax_fit.c
trunk/target_functions/relax_fit.h
Modified: trunk/target_functions/exponential.c
URL:
http://svn.gna.org/viewcvs/relax/trunk/target_functions/exponential.c?rev=25249&r1=25248&r2=25249&view=diff
==============================================================================
--- trunk/target_functions/exponential.c (original)
+++ trunk/target_functions/exponential.c Mon Aug 25 15:28:43 2014
@@ -21,6 +21,9 @@
/* The exponential function is needed. */
#include <math.h>
+/* functions for the exponential */
+#include "exponential.h"
+
void exponential(double *params, double *relax_times, double *back_calc, int
num_times) {
/* Function to back calculate the peak intensities.
@@ -43,3 +46,46 @@
}
}
+
+void exponential_dI(double *params, double *relax_times, double
back_calc_grad[][MAXTIMES], int num_times) {
+ /* Calculate the dI partial derivate of the 2-parameter exponential
curve.
+ */
+
+ /* Declarations */
+ 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_grad[1][i] = 0.0;
+
+ /* The partial derivate */
+ else
+ back_calc_grad[1][i] = exp(-relax_times[i] * params[0]);
+ }
+}
+
+
+void exponential_dR(double *params, double *relax_times, double
back_calc_grad[][MAXTIMES], int num_times) {
+ /* Calculate the dR partial derivate of the 2-parameter exponential
curve.
+ */
+
+ /* Declarations */
+ 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_grad[0][i] = 0.0;
+
+ /* The partial derivate */
+ else
+ back_calc_grad[0][i] = -params[1] * relax_times[i] *
exp(-relax_times[i] * params[0]);
+ }
+}
Modified: trunk/target_functions/exponential.h
URL:
http://svn.gna.org/viewcvs/relax/trunk/target_functions/exponential.h?rev=25249&r1=25248&r2=25249&view=diff
==============================================================================
--- trunk/target_functions/exponential.h (original)
+++ trunk/target_functions/exponential.h Mon Aug 25 15:28:43 2014
@@ -1,6 +1,7 @@
/*
* Copyright (C) 2006 Gary S Thompson (see https://gna.org/users for contact
* details)
+ * Copyright (C) 2014 Edward d'Auvergne
*
* This file is part of the program relax (http://www.nmr-relax.com).
*
@@ -20,6 +21,12 @@
#ifndef RELAX_EXPONENTIAL
#define RELAX_EXPONENTIAL
+/* The maximum number of spectral time points */
+#define MAXTIMES 50
+
+
void exponential(double *params, double *relax_times, double *back_calc, int
num_times);
+void exponential_dI(double *params, double *relax_times, double
back_calc_grad[][MAXTIMES], int num_times);
+void exponential_dR(double *params, double *relax_times, double
back_calc_grad[][MAXTIMES], int num_times);
#endif
Modified: trunk/target_functions/relax_fit.c
URL:
http://svn.gna.org/viewcvs/relax/trunk/target_functions/relax_fit.c?rev=25249&r1=25248&r2=25249&view=diff
==============================================================================
--- trunk/target_functions/relax_fit.c (original)
+++ trunk/target_functions/relax_fit.c Mon Aug 25 15:28:43 2014
@@ -1,5 +1,5 @@
/*
- * Copyright (C) 2006-2013 Edward d'Auvergne
+ * Copyright (C) 2006-2014 Edward d'Auvergne
*
* This file is part of the program relax (http://www.nmr-relax.com).
*
@@ -164,6 +164,54 @@
/* Return the numpy array */
return back_calc_py;
+}
+
+
+static PyObject *
+jacobian(PyObject *self, PyObject *args) {
+ /* Return the Jacobian as a Python list. */
+
+ /* Declarations */
+ PyObject *params_arg;
+ PyObject *element;
+ int i, j;
+
+ /* The Python list of lists. */
+ PyObject *list = PyList_New(num_params);
+ Py_INCREF(list);
+
+ /* Parse the function arguments, the only argument should be the
parameter array */
+ if (!PyArg_ParseTuple(args, "O", ¶ms_arg))
+ return NULL;
+
+ /* Place the parameter array elements into the C array */
+ for (i = 0; i < num_params; i++) {
+ /* Get the element */
+ element = PySequence_GetItem(params_arg, i);
+
+ /* Convert to a C double, then free the memory. */
+ params[i] = PyFloat_AsDouble(element);
+ Py_CLEAR(element);
+
+ /* Scale the parameter */
+ params[i] = params[i] * scaling_matrix[i];
+ }
+
+ /* The partial derivates */
+ exponential_dR(params, relax_times, back_calc_grad, num_times);
+ exponential_dI(params, relax_times, back_calc_grad, num_times);
+
+ /* Convert to a Python list */
+ for (i = 0; i < num_params; i++) {
+ PyObject *list2 = PyList_New(num_times);
+ for (j = 0; j < num_times; j++) {
+ PyList_Append(list2, PyFloat_FromDouble(back_calc_grad[i][j]));
+ }
+ PyList_Append(list, list2);
+ }
+
+ /* Return the Jacobian */
+ return list;
}
@@ -194,6 +242,11 @@
back_calc_I,
METH_VARARGS,
"Return the back calculated peak intensities as a Python list."
+ }, {
+ "jacobian",
+ jacobian,
+ METH_VARARGS,
+ "Return the Jacobian matrix as a Python list."
},
{NULL, NULL, 0, NULL} /* Sentinel */
};
Modified: trunk/target_functions/relax_fit.h
URL:
http://svn.gna.org/viewcvs/relax/trunk/target_functions/relax_fit.h?rev=25249&r1=25248&r2=25249&view=diff
==============================================================================
--- trunk/target_functions/relax_fit.h (original)
+++ trunk/target_functions/relax_fit.h Mon Aug 25 15:28:43 2014
@@ -39,6 +39,7 @@
/* Variables used for storage during the function calls of optimisation */
static double back_calc[MAXTIMES];
+static double back_calc_grad[MAXPARAMS][MAXTIMES];
static double params[MAXPARAMS];
static double values[MAXTIMES];
static double sd[MAXTIMES];
r25249 - in /trunk/target_functions: exponential.c exponential.h relax_fit.c relax_fit.h