r25303 - in /trunk/target_functions: c_chi2.c c_chi2.h -- August 26, 2014 - 18:08

Author: bugman
Date: Tue Aug 26 18:08:15 2014
New Revision: 25303

URL: http://svn.gna.org/viewcvs/relax?rev=25303&view=rev
Log:
Implemented the C version of the chi-squared Hessian.

This is a direct translation of the Python code.


Modified:
    trunk/target_functions/c_chi2.c
    trunk/target_functions/c_chi2.h

Modified: trunk/target_functions/c_chi2.c
URL: 
http://svn.gna.org/viewcvs/relax/trunk/target_functions/c_chi2.c?rev=25303&r1=25302&r2=25303&view=diff
==============================================================================
--- trunk/target_functions/c_chi2.c     (original)
+++ trunk/target_functions/c_chi2.c     Tue Aug 26 18:08:15 2014
@@ -113,3 +113,63 @@
         }
     }
 }
+
+
+void d2chi2(double d2chi2[MAX_PARAMS][MAX_PARAMS], double data[MAX_DATA], 
double back_calc_vals[MAX_DATA], double back_calc_grad[MAX_PARAMS][MAX_DATA], 
double back_calc_hess[MAX_PARAMS][MAX_PARAMS][MAX_DATA], double 
errors[MAX_DATA], int num_points, int num_params) {
+    /* Calculate the full chi-squared Hessian.
+
+    The chi-squared Hessian
+    =======================
+
+    The equation is::
+
+                              _n_
+        d2chi^2(theta)        \       1      / dyi(theta)   dyi(theta)       
                 d2yi(theta)   \
+        ---------------  =  2  >  ---------- | ---------- . ----------  -  
(yi-yi(theta)) . --------------- |
+        dthetaj.dthetak       /__ sigma_i**2 \  dthetaj      dthetak         
               dthetaj.dthetak /
+                              i=1
+
+    where
+        - i is the index over data sets.
+        - j is the parameter index for the first dimension of the Hessian.
+        - k is the parameter index for the second dimension of the Hessian.
+        - 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.
+        - d2yi(theta)/dthetaj.dthetak are the values of the back calculated 
Hessian for the
+        parameters j and k.
+        - sigma_i are the values of the error set.
+
+
+    @param d2chi2:              The chi-squared Hessian data structure to 
place the Hessian elements into.
+    @type d2chi2:               numpy rank-2 size MxM array
+    @param data:                The vector of yi values.
+    @type data:                 numpy rank-1 size N array
+    @param back_calc_vals:      The vector of yi(theta) values.
+    @type back_calc_vals:       numpy rank-1 size N array
+    @param back_calc_grad:      The matrix of dyi(theta)/dtheta values.
+    @type back_calc_grad:       numpy rank-2 size MxN array
+    @param back_calc_hess:      The matrix of d2yi(theta)/dtheta.dtheta 
values.
+    @type back_calc_hess:       numpy rank-3 size MxMxN array
+    @param errors:              The vector of sigma_i values.
+    @type errors:               numpy rank-1 size N array
+    @param num_points:          The number of data points to sum over.
+    @type num_points:           int
+    @param num_params:          The dimensions of the Hessian.
+    @type num_params:           int
+    */
+
+    /* Declarations. */
+    int i, j, k;
+
+    /* Calculate the chi-squared Hessian. */
+    for (j = 0; j < num_params; ++j) {
+        for (k = 0; k < num_params; ++k) {
+            d2chi2[j][k] = 0.0;
+            for (i = 0; i < num_points; ++i) {
+                d2chi2[j][k] += 2.0 / square(errors[i]) * 
(back_calc_grad[j][i] * back_calc_grad[k][i] - (data[i] - back_calc_vals[i]) 
* back_calc_hess[j][k][i]);
+            }
+        }
+    }
+}

Modified: trunk/target_functions/c_chi2.h
URL: 
http://svn.gna.org/viewcvs/relax/trunk/target_functions/c_chi2.h?rev=25303&r1=25302&r2=25303&view=diff
==============================================================================
--- trunk/target_functions/c_chi2.h     (original)
+++ trunk/target_functions/c_chi2.h     Tue Aug 26 18:08:15 2014
@@ -22,10 +22,12 @@
 #ifndef RELAX_C_CHI2 
 #define RELAX_C_CHI2
 
-/* The maximum number of data points */
+/* The maximum number of parameters and data points */
+#define MAX_PARAMS 10
 #define MAX_DATA 50
 
 double chi2(double values[], double sd[], double back_calc[], int num_times);
 void dchi2(double dchi2[], double data[], double back_calc_vals[], double 
back_calc_grad[][MAX_DATA], double errors[], int num_times, int M);
+void d2chi2(double d2chi2[MAX_PARAMS][MAX_PARAMS], double data[MAX_DATA], 
double back_calc_vals[MAX_DATA], double back_calc_grad[MAX_PARAMS][MAX_DATA], 
double back_calc_hess[MAX_PARAMS][MAX_PARAMS][MAX_DATA], double 
errors[MAX_DATA], int num_times, int M);
 
 #endif