diff options
Diffstat (limited to 'src/post-refinement.c')
-rw-r--r-- | src/post-refinement.c | 334 |
1 files changed, 6 insertions, 328 deletions
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<num_params; k++ ) { - gradients[k] = gradient(cr, rv[k], refl, pmodel); - } - - for ( k=0; k<num_params; k++ ) { - - int g; - double v_c, v_curr; - - for ( g=0; g<num_params; g++ ) { - - double M_c, M_curr; - - /* Matrix is symmetric */ - if ( g > 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<num_params; param++ ) { - double shift = gsl_vector_get(shifts, param); - apply_shift(cr, rv[param], shift); - if ( fabs(shift) > 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<n_crystals; i++ ) { - - double r, free_r, log_r, free_log_r; - - if ( crystal_get_user_flag(crystals[i]) ) continue; - - r = residual(crystals[i], full, 0, NULL, NULL); - free_r = residual(crystals[i], full, 1, NULL, NULL); - log_r = log_residual(crystals[i], full, 0, NULL, NULL); - free_log_r = log_residual(crystals[i], full, 1, NULL, NULL); - - if ( isnan(r) || isnan(free_r) - || isnan(log_r) || isnan(free_log_r) ) continue; - - *presidual += r; - *pfree_residual += free_r; - *plog_residual += log_r; - *pfree_log_residual += free_log_r; - } -} - - static void write_residual_graph(Crystal *cr, const RefList *full) { int i; @@ -877,49 +619,6 @@ static void write_residual_graph(Crystal *cr, const RefList *full) } -static void do_scale_refine(Crystal *cr, const RefList *full, - PartialityModel pmodel, int verbose) -{ - int i, done; - double old_dev; - - old_dev = log_residual(cr, full, 0, NULL, NULL); - - if ( verbose ) { - STATUS("Initial G=%.2f, B=%e\n", - crystal_get_osf(cr), crystal_get_Bfac(cr)); - STATUS("Scaling initial dev = %10.5e, free dev = %10.5e\n", - old_dev, log_residual(cr, full, 1, NULL, NULL)); - } - - i = 0; - done = 0; - do { - - double dev; - - scale_iterate(cr, full, pmodel); - - dev = log_residual(cr, full, 0, 0, NULL); - if ( fabs(dev - old_dev) < dev*0.01 ) done = 1; - - if ( verbose ) { - STATUS("Scaling %2i: dev = %10.5e, free dev = %10.5e\n", - i+1, dev, log_residual(cr, full, 0, 0, NULL)); - } - - i++; - old_dev = dev; - - } while ( i < MAX_CYCLES && !done ); - - if ( verbose ) { - STATUS("Final G=%.2f, B=%e\n", crystal_get_osf(cr), - crystal_get_Bfac(cr)); - } -} - - static void do_pr_refine(Crystal *cr, const RefList *full, PartialityModel pmodel, int verbose) { @@ -977,8 +676,7 @@ static void do_pr_refine(Crystal *cr, const RefList *full, static struct prdata pr_refine(Crystal *cr, const RefList *full, - PartialityModel pmodel, int no_scale, int no_pr, - double max_B) + PartialityModel pmodel) { int verbose = 0; struct prdata prdata; @@ -986,23 +684,11 @@ static struct prdata pr_refine(Crystal *cr, const RefList *full, prdata.refined = 0; prdata.n_filtered = 0; - if ( !no_scale ) { - do_scale_refine(cr, full, pmodel, verbose); - } - - /* Reject if B factor modulus is very large */ - if ( fabs(crystal_get_Bfac(cr)) > 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; |