Author: bugman
Date: Fri Aug 29 09:06:01 2014
New Revision: 25413
URL: http://svn.gna.org/viewcvs/relax?rev=25413&view=rev
Log:
First attempt at properly implementing the
target_functions.relax_fit.jacobian() function.
This is now the Jacobian of the chi-squared function. A new jacobian_matrix
data structure has been
created for holding the matrix data prior to converting it into a Python list
of lists. The
equation used was simply the chi-squared gradient whereby the sum over i has
been dropped and the i
elements are stored in the second dimension of matrix.
Modified:
trunk/target_functions/relax_fit.c
trunk/target_functions/relax_fit.h
Modified: trunk/target_functions/relax_fit.c
URL:
http://svn.gna.org/viewcvs/relax/trunk/target_functions/relax_fit.c?rev=25413&r1=25412&r2=25413&view=diff
==============================================================================
--- trunk/target_functions/relax_fit.c (original)
+++ trunk/target_functions/relax_fit.c Fri Aug 29 09:06:01 2014
@@ -232,7 +232,27 @@
static PyObject *
jacobian(PyObject *self, PyObject *args) {
- /* Return the Jacobian as a Python list of lists. */
+ /* Return the Jacobian as a Python list of lists.
+
+ The Jacobian
+ ============
+
+ The equation is::
+
+ / yi - yi(theta) dyi(theta) \
+ J_ji = -2 | -------------- . ---------- |
+ \ sigma_i**2 dthetaj /
+
+ where
+ - i is the index over data sets.
+ - j is the parameter index.
+ - theta is the parameter vector.
+ - yi are the values of the measured data set.
+ - yi(theta) are the values of the back calculated data set.
+ - dyi(theta)/dthetaj are the values of the back calculated gradient
for parameter j.
+ - sigma_i are the values of the error set.
+
+ */
/* Declarations. */
PyObject *params_arg;
@@ -246,9 +266,19 @@
/* Convert the parameters Python list to a C array. */
param_to_c(params_arg);
+ /* Back calculated the peak intensities. */
+ exponential(params[index_I0], params[index_R], relax_times, back_calc,
num_times);
+
/* The partial derivatives. */
exponential_dR(params[index_I0], params[index_R], index_R, relax_times,
back_calc_grad, num_times);
exponential_dI0(params[index_I0], params[index_R], index_I0, relax_times,
back_calc_grad, num_times);
+
+ /* Assemble the chi-squared Jacobian. */
+ for (j = 0; j < num_params; ++j) {
+ for (i = 0; i < num_times; ++i) {
+ jacobian_matrix[j][i] = -2.0 / square(sd[i]) * (values[i] -
back_calc[i]) * back_calc_grad[j][i];
+ }
+ }
/* Convert to a Python list of lists. */
list = PyList_New(0);
@@ -257,7 +287,7 @@
list2 = PyList_New(0);
Py_INCREF(list2);
for (j = 0; j < num_times; j++) {
- PyList_Append(list2, PyFloat_FromDouble(back_calc_grad[i][j]));
+ PyList_Append(list2, PyFloat_FromDouble(jacobian_matrix[i][j]));
}
PyList_Append(list, list2);
}
Modified: trunk/target_functions/relax_fit.h
URL:
http://svn.gna.org/viewcvs/relax/trunk/target_functions/relax_fit.h?rev=25413&r1=25412&r2=25413&view=diff
==============================================================================
--- trunk/target_functions/relax_fit.h (original)
+++ trunk/target_functions/relax_fit.h Fri Aug 29 09:06:01 2014
@@ -26,6 +26,9 @@
#define PyMODINIT_FUNC void
#endif
+/* Define the function for calculating the square of a number. */
+#define square(x) ((x)*(x))
+
/****************************************/
/* External, hence permanent, variables. */
@@ -44,6 +47,7 @@
static double back_calc_hess[MAX_PARAMS][MAX_PARAMS][MAX_DATA];
static double dchi2_vals[MAX_PARAMS];
static double d2chi2_vals[MAX_PARAMS][MAX_PARAMS];
+static double jacobian_matrix[MAX_PARAMS][MAX_DATA];
static double params[MAX_PARAMS];
static double values[MAX_DATA];
static double sd[MAX_DATA];
_______________________________________________
relax (http://www.nmr-relax.com)
This is the relax-commits mailing list
[email protected]
To unsubscribe from this list, get a password
reminder, or change your subscription options,
visit the list information page at
https://mail.gna.org/listinfo/relax-commits
