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
_______________________________________________
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
