aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2018-01-31 15:39:17 +0100
committerThomas White <taw@physics.org>2018-02-27 17:12:42 +0100
commite0066f57d2d9b4b86d7d1919f8adc71c2e6a2bdb (patch)
treece3ab9b8a006c7c1dbebfb9351aaaae026a5fa31
parenta803b5205fb9fddf2d94036999522394733a3732 (diff)
Stop if GSL gives us NAN residual
Happens if the initial position is close to a tight minimum, so all simplex vertices have bad values.
-rw-r--r--src/post-refinement.c83
1 files changed, 54 insertions, 29 deletions
diff --git a/src/post-refinement.c b/src/post-refinement.c
index ca550782..088e25dd 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -457,11 +457,19 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
}
do {
+ double res;
+
n_iter++;
status = gsl_multimin_fminimizer_iterate(min);
if ( status ) break;
+ res = residual_f(min->x, &residual_f_priv);
+ if ( isnan(res) ) {
+ status = GSL_ENOPROG;
+ break;
+ }
+
if ( verbose ) {
double res = residual_f(min->x, &residual_f_priv);
double size = gsl_multimin_fminimizer_size(min);
@@ -495,42 +503,59 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
STATUS("status = %i (%s)\n", status, gsl_strerror(status));
}
- if ( check_angle_shifts(min->x, initial, rv, n_params, &residual_f_priv) ) return;
+ if ( status == GSL_SUCCESS ) {
- if ( verbose ) {
+ if ( check_angle_shifts(min->x, initial, rv, n_params, &residual_f_priv) ) return;
- double res = residual_f(min->x, &residual_f_priv);
- double size = gsl_multimin_fminimizer_size(min);
- STATUS("At end: %f %f %f %f ----> %f %f %e %f residual = %e size %f\n",
- gsl_vector_get(min->x, 0),
- gsl_vector_get(min->x, 1),
- gsl_vector_get(min->x, 2),
- gsl_vector_get(min->x, 3),
- rad2deg(get_actual_val(min->x, initial, rv, 0)),
- rad2deg(get_actual_val(min->x, initial, rv, 1)),
- get_actual_val(min->x, initial, rv, 2),
- get_actual_val(min->x, initial, rv, 3)*1e10,
- res, size);
+ if ( verbose ) {
- }
+ double res = residual_f(min->x, &residual_f_priv);
+ double size = gsl_multimin_fminimizer_size(min);
+ STATUS("At end: %f %f %f %f ----> %f %f %e %f residual = %e size %f\n",
+ gsl_vector_get(min->x, 0),
+ gsl_vector_get(min->x, 1),
+ gsl_vector_get(min->x, 2),
+ gsl_vector_get(min->x, 3),
+ rad2deg(get_actual_val(min->x, initial, rv, 0)),
+ rad2deg(get_actual_val(min->x, initial, rv, 1)),
+ get_actual_val(min->x, initial, rv, 2),
+ get_actual_val(min->x, initial, rv, 3)*1e10,
+ res, size);
- /* Apply the final shifts */
- apply_parameters(min->x, initial, rv, cr);
- update_predictions(cr);
- calculate_partialities(cr, PMODEL_XSPHERE);
- r = linear_scale(full, crystal_get_reflections(cr), &G, 0);
- if ( r == 0 ) {
- crystal_set_osf(cr, G);
- }
+ }
- if ( verbose ) {
+ if ( fh != NULL ) {
+ double res = residual_f(min->x, &residual_f_priv);
+ fprintf(fh, "%5i %10.8f %10.8f %5i %10.8f %10.8f %e %e\n",
+ n_iter, res, 0.0, 0,
+ rad2deg(get_actual_val(min->x, initial, rv, 0)),
+ rad2deg(get_actual_val(min->x, initial, rv, 1)),
+ get_actual_val(min->x, initial, rv, 2),
+ get_actual_val(min->x, initial, rv, 3)*1e10);
+ fclose(fh);
+ }
- STATUS("After applying final shifts:\n");
- STATUS("PR final: dev = %10.5e, free dev = %10.5e\n",
- residual(cr, full, 0, NULL, NULL, 0),
- residual(cr, full, 1, NULL, NULL, 0));
- STATUS("Final R = %e m^-1\n", crystal_get_profile_radius(cr));
+ /* Apply the final shifts */
+ apply_parameters(min->x, initial, rv, cr);
+ update_predictions(cr);
+ calculate_partialities(cr, PMODEL_XSPHERE);
+ r = linear_scale(full, crystal_get_reflections(cr), &G, 0);
+ if ( r == 0 ) {
+ crystal_set_osf(cr, G);
+ }
+
+ if ( verbose ) {
+
+ STATUS("After applying final shifts:\n");
+ STATUS("PR final: dev = %10.5e, free dev = %10.5e\n",
+ residual(cr, full, 0, NULL, NULL, 0),
+ residual(cr, full, 1, NULL, NULL, 0));
+ STATUS("Final R = %e m^-1\n", crystal_get_profile_radius(cr));
+
+ }
+ } else {
+ ERROR("Bad refinement: crystal %i\n", serial);
}
gsl_multimin_fminimizer_free(min);