Some code simplification

hrue · Dec 20, 2023 · e9d07e8 · e9d07e8
1 parent 0a85b8c
commit e9d07e8
Show file tree

Hide file tree

Showing 2 changed files with 44 additions and 24 deletions.
diff --git a/gmrflib/blockupdate.c b/gmrflib/blockupdate.c
@@ -213,27 +213,50 @@ int GMRFLib_2order_approx_core(int thread_id, double *a, double *b, double *c, d
 		df = (1.0 / 12.0 * f[4] - 2.0 / 3.0 * f[3] + 0.0 * f[2] + 2.0 / 3.0 * f[1] - 1.0 / 12.0 * f[0]) / step;
 		ddf = (-1.0 / 12.0 * f[4] - 4.0 / 3.0 * f[3] - 5.0 / 2.0 * f[2] + 4.0 / 3.0 * f[1] - 1.0 / 12.0 * f[0]) / SQR(step);
 		dddf = (-1.0 / 2.0 * f[4] + 1.0 * f[3] + 0.0 * f[2] - 1.0 * f[1] + 1.0 / 2.0 * f[0]) / POW3(step);
+	} else if (stencil && *stencil == 3) {
+		// special implementation for *stencil==3, which is used for the Gaussian case
+		if (!step_len || ISZERO(*step_len)) {
+			step = 1.0e-4 * (GSL_DBL_EPSILON / 2.220446049e-16);
+		} else {
+			step = *step_len;
+		}
+
+		int n = 3;
+		xx[0] = x0 - step;
+		xx[1] = x0;
+		xx[2] = x0 + step;
+
+		loglFunc(thread_id, f, xx, n, indx, x_vec, NULL, loglFunc_arg, NULL);
+
+		f0 = f[1];
+		df = 0.5 * (- f[0] + f[2]);
+		ddf = f[0] - 2.0 * f[1] + f[2];
+
+		df /= step;
+		ddf /= SQR(step);
+		dddf = 0.0;
 	} else {
 		int num_points = (stencil ? *stencil : 5);
 
 		if (!step_len || ISZERO(*step_len)) {
-			step = (GSL_DBL_EPSILON / 2.220446049e-16) * (*stencil == 5 ? 1.0e-4 : (*stencil == 7 ? 5.0e-4 : 1.0e-3));
+			step = (GSL_DBL_EPSILON / 2.220446049e-16) * (*stencil == 3 || *stencil == 5 ? 1.0e-4 : (*stencil == 7 ? 5.0e-4 : 1.0e-3));
 		} else {
 			step = *step_len;
 		}
 
 		// see https://en.wikipedia.org/wiki/Finite_difference_coefficients
 		int n, nn, wlength;
 
-		static double wf3[] = {
-			// 
-			-0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0,
-			// 
-			1.0, -2.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
-			// 
-			0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
-		};
-
+		/*
+		  static double wf3[] = {
+		  // 
+		  -0.5, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0, 0.0,
+		  // 
+		  1.0, -2.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
+		  // 
+		  0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
+		  };
+		*/
 		static double wf5[] = {
 			// 
 			0.08333333333333333, -0.6666666666666666, 0.0, 0.6666666666666666, -0.08333333333333333, 0.0, 0.0, 0.0,
@@ -267,12 +290,14 @@ int GMRFLib_2order_approx_core(int thread_id, double *a, double *b, double *c, d
 		double *wf = NULL;
 
 		switch (num_points) {
-		case 3:
-			n = 3;
-			nn = 2;
-			wlength = 8;
-			wf = wf3;
-			break;
+			/*
+			  case 3:
+			  n = 3;
+			  nn = 1;
+			  wlength = 8;
+			  wf = wf3;
+			  break;
+			*/
 		case 5:
 			n = 5;
 			nn = 2;

diff --git a/inlaprog/src/inla-likelihood.c b/inlaprog/src/inla-likelihood.c
@@ -663,18 +663,14 @@ int loglikelihood_gaussian(int thread_id, double *logll, double *x, int m, int i
 	}
 	Data_section_tp *ds = (Data_section_tp *) arg;
 	double y, lprec, prec, w;
-	static double log_prec_limit = 0.0;
+	static double log_prec_limit = -log(INLA_REAL_SMALL);
+
 	y = ds->data_observations.y[idx];
 	w = ds->data_observations.weight_gaussian[idx];
 
-	if (log_prec_limit == 0.0) {
-		log_prec_limit = -log(INLA_REAL_SMALL);
-	}
-
 	if (ds->data_observations.log_prec_gaussian_offset[thread_id][0] > log_prec_limit) {
 		lprec = ds->data_observations.log_prec_gaussian[thread_id][0] + log(w);
 		prec = exp(lprec);
-		// prec = map_precision(ds->data_observations.log_prec_gaussian[thread_id][0], MAP_FORWARD, NULL) * w;
 	} else {
 		double prec_offset = map_precision(ds->data_observations.log_prec_gaussian_offset[thread_id][0], MAP_FORWARD, NULL);
 		double prec_var = map_precision(ds->data_observations.log_prec_gaussian[thread_id][0], MAP_FORWARD, NULL);
@@ -697,8 +693,7 @@ int loglikelihood_gaussian(int thread_id, double *logll, double *x, int m, int i
 
 	if (m > 0) {
 		if (PREDICTOR_LINK_EQ(link_identity)) {
-
-			if (PREDICTOR_LINK_EQ(link_identity) && (PREDICTOR_SCALE == 1.0 && off == 0.0)) {
+			if (PREDICTOR_SCALE == 1.0 && off == 0.0) {
 				double a = -0.5 * prec;
 				double b = LOG_NORMC_GAUSSIAN + 0.5 * lprec;
 				if (0 && m >= 8L) {