From a2b11362e440d3487520a39284ae72fcb4cb37f5 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 11 Nov 2015 17:31:34 +0100 Subject: partialator: Scale (strictly) using strong reflections only --- src/post-refinement.c | 334 +------------------------------------------------- 1 file changed, 6 insertions(+), 328 deletions(-) (limited to 'src/post-refinement.c') diff --git a/src/post-refinement.c b/src/post-refinement.c index cf3108db..487ce872 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -216,21 +216,9 @@ double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) struct image *image = crystal_get_image(cr); double R = crystal_get_profile_radius(cr); double gr; - signed int hi, ki, li; - double s; - get_indices(refl, &hi, &ki, &li); - s = resolution(crystal_get_cell(cr), hi, ki, li); get_partial(refl, &rlow, &rhigh, &p); - if ( k == GPARAM_OSF ) { - return 1.0; - } - - if ( k == GPARAM_BFAC ) { - return -s*s; - } - if ( k == GPARAM_R ) { double Rglow, Rghigh; @@ -347,18 +335,6 @@ static void apply_shift(Crystal *cr, int k, double shift) crystal_set_profile_radius(cr, t); break; - case GPARAM_BFAC : - t = crystal_get_Bfac(cr); - t += shift; - crystal_set_Bfac(cr, t); - break; - - case GPARAM_OSF : - t = crystal_get_osf(cr); - t += shift; - crystal_set_osf(cr, t); - break; - case GPARAM_ASX : case GPARAM_ASY : case GPARAM_ASZ : @@ -379,141 +355,6 @@ static void apply_shift(Crystal *cr, int k, double shift) } -/* Perform one cycle of scaling of 'cr' against 'full' */ -static double scale_iterate(Crystal *cr, const RefList *full, - PartialityModel pmodel) -{ - gsl_matrix *M; - gsl_vector *v; - gsl_vector *shifts; - int param; - Reflection *refl; - RefListIterator *iter; - RefList *reflections; - double max_shift; - int nref = 0; - int num_params = 0; - enum gparam rv[32]; - double G, B; - - rv[num_params++] = GPARAM_OSF; - rv[num_params++] = GPARAM_BFAC; - - M = gsl_matrix_calloc(num_params, num_params); - v = gsl_vector_calloc(num_params); - - reflections = crystal_get_reflections(cr); - G = crystal_get_osf(cr); - B = crystal_get_Bfac(cr); - - /* Scaling terms */ - for ( refl = first_refl(reflections, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - signed int ha, ka, la; - double I_full, delta_I, esd; - double w; - double I_partial; - int k; - double p, L, s; - double fx; - Reflection *match; - double gradients[num_params]; - - /* If reflection is free-flagged, don't use it here */ - if ( get_flag(refl) ) continue; - - /* Find the full version */ - get_indices(refl, &ha, &ka, &la); - match = find_refl(full, ha, ka, la); - if ( match == NULL ) continue; - - /* Merged intensitty */ - I_full = get_intensity(match); - - /* Actual measurement of this reflection from this pattern */ - I_partial = get_intensity(refl); - esd = get_esd_intensity(refl); - p = get_partiality(refl); - - /* Scale only using strong reflections */ - if ( I_partial <= 3.0*esd ) continue; /* Also because of log */ - if ( get_redundancy(match) < 2 ) continue; - if ( I_full <= 0 ) continue; /* Because log */ - if ( p <= 0.0 ) continue; /* Because of log */ - - L = get_lorentz(refl); - s = resolution(crystal_get_cell(cr), ha, ka, la); - - /* Calculate the weight for this reflection */ - w = 1.0; - - /* Calculate all gradients for this reflection */ - for ( k=0; k k ) continue; - - M_c = w * gradients[g] * gradients[k]; - - M_curr = gsl_matrix_get(M, k, g); - gsl_matrix_set(M, k, g, M_curr + M_c); - gsl_matrix_set(M, g, k, M_curr + M_c); - - } - - fx = G + log(p) - log(L) - B*s*s + log(I_full); - delta_I = log(I_partial) - fx; - v_c = w * delta_I * gradients[k]; - v_curr = gsl_vector_get(v, k); - gsl_vector_set(v, k, v_curr + v_c); - - } - - nref++; - } - - if ( nref < num_params ) { - crystal_set_user_flag(cr, PRFLAG_FEWREFL); - gsl_matrix_free(M); - gsl_vector_free(v); - return 0.0; - } - - max_shift = 0.0; - shifts = solve_svd(v, M, NULL, 0); - if ( shifts != NULL ) { - - for ( param=0; param max_shift ) max_shift = fabs(shift); - } - - } else { - crystal_set_user_flag(cr, PRFLAG_SOLVEFAIL); - } - - gsl_matrix_free(M); - gsl_vector_free(v); - gsl_vector_free(shifts); - - return max_shift; -} - - /* Perform one cycle of post refinement on 'image' against 'full' */ static double pr_iterate(Crystal *cr, const RefList *full, PartialityModel pmodel, @@ -667,74 +508,8 @@ static double pr_iterate(Crystal *cr, const RefList *full, } -static double log_residual(Crystal *cr, const RefList *full, int free, - int *pn_used, const char *filename) -{ - double dev = 0.0; - double G, B; - Reflection *refl; - RefListIterator *iter; - int n_used = 0; - FILE *fh = NULL; - - G = crystal_get_osf(cr); - B = crystal_get_Bfac(cr); - if ( filename != NULL ) { - fh = fopen(filename, "a"); - if ( fh == NULL ) { - ERROR("Failed to open '%s'\n", filename); - } - } - - for ( refl = first_refl(crystal_get_reflections(cr), &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - double p, L, s, w; - signed int h, k, l; - Reflection *match; - double esd, I_full, I_partial; - double fx, dc; - - if ( free && !get_flag(refl) ) continue; - - get_indices(refl, &h, &k, &l); - match = find_refl(full, h, k, l); - if ( match == NULL ) continue; - - p = get_partiality(refl); - L = get_lorentz(refl); - I_partial = get_intensity(refl); - I_full = get_intensity(match); - esd = get_esd_intensity(refl); - s = resolution(crystal_get_cell(cr), h, k, l); - - if ( I_partial <= 3.0*esd ) continue; /* Also because of log */ - if ( get_redundancy(match) < 2 ) continue; - if ( I_full <= 0 ) continue; /* Because log */ - if ( p <= 0.0 ) continue; /* Because of log */ - - fx = G + log(p) - log(L) - B*s*s + log(I_full); - dc = log(I_partial) - fx; - w = 1.0; - dev += w*dc*dc; - - if ( fh != NULL ) { - fprintf(fh, "%4i %4i %4i %e %e\n", - h, k, l, s, dev); - } - - } - - if ( fh != NULL ) fclose(fh); - - if ( pn_used != NULL ) *pn_used = n_used; - return dev; -} - - -static double residual(Crystal *cr, const RefList *full, int free, - int *pn_used, const char *filename) +double residual(Crystal *cr, const RefList *full, int free, + int *pn_used, const char *filename) { double dev = 0.0; double G, B; @@ -799,39 +574,6 @@ static double residual(Crystal *cr, const RefList *full, int free, } -void all_residuals(Crystal **crystals, int n_crystals, RefList *full, - double *presidual, double *pfree_residual, - double *plog_residual, double *pfree_log_residual) -{ - int i; - - *presidual = 0.0; - *pfree_residual = 0.0; - *plog_residual = 0.0; - *pfree_log_residual = 0.0; - - for ( i=0; i max_B ) { - crystal_set_user_flag(cr, PRFLAG_BIGB); - return prdata; - } - if ( verbose ) { write_residual_graph(cr, full); } - if ( !no_pr ) { - do_pr_refine(cr, full, pmodel, verbose); - } + do_pr_refine(cr, full, pmodel, verbose); if ( crystal_get_user_flag(cr) == 0 ) { prdata.refined = 1; @@ -1017,9 +703,6 @@ struct refine_args RefList *full; Crystal *crystal; PartialityModel pmodel; - int no_scale; - int no_pr; - double max_B; struct prdata prdata; }; @@ -1039,8 +722,7 @@ static void refine_image(void *task, int id) struct refine_args *pargs = task; Crystal *cr = pargs->crystal; - pargs->prdata = pr_refine(cr, pargs->full, pargs->pmodel, - pargs->no_scale, pargs->no_pr, pargs->max_B); + pargs->prdata = pr_refine(cr, pargs->full, pargs->pmodel); } @@ -1072,8 +754,7 @@ static void done_image(void *vqargs, void *task) void refine_all(Crystal **crystals, int n_crystals, - RefList *full, int nthreads, PartialityModel pmodel, - int no_scale, int no_pr, double max_B) + RefList *full, int nthreads, PartialityModel pmodel) { struct refine_args task_defaults; struct queue_args qargs; @@ -1083,9 +764,6 @@ void refine_all(Crystal **crystals, int n_crystals, task_defaults.pmodel = pmodel; task_defaults.prdata.refined = 0; task_defaults.prdata.n_filtered = 0; - task_defaults.no_scale = no_scale; - task_defaults.no_pr = no_pr; - task_defaults.max_B = max_B; qargs.task_defaults = task_defaults; qargs.n_started = 0; -- cgit v1.2.3