r25413 - in /trunk/target_functions: relax_fit.c relax_fit.h -- August 29, 2014 - 09:06

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];