diff options
| -rw-r--r-- | src/partialator.c | 13 | ||||
| -rw-r--r-- | src/post-refinement.c | 74 | ||||
| -rw-r--r-- | src/post-refinement.h | 18 |
3 files changed, 78 insertions, 27 deletions
diff --git a/src/partialator.c b/src/partialator.c index 079d4d46..5032296f 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -991,19 +991,28 @@ int main(int argc, char *argv[]) for ( i=0; i<n_iter; i++ ) { double init_dev, init_free_dev; + double init_log_dev, init_free_log_dev; double final_dev, final_free_dev; + double final_log_dev, final_free_log_dev; STATUS("Refinement cycle %i of %i\n", i+1, n_iter); /* Refine all crystals to get the best fit */ refine_all(crystals, n_crystals, full, nthreads, pmodel, - no_scale, no_pr, &init_dev, &init_free_dev, - &final_dev, &final_free_dev); + no_scale, no_pr, + &init_dev, &init_free_dev, + &init_log_dev, &init_free_log_dev, + &final_dev, &final_free_dev, + &final_log_dev, &final_free_log_dev); STATUS("Overall residual: initial = %e, final = %e\n", init_dev, final_dev); STATUS("Overall free residual: initial = %e, final = %e\n", init_free_dev, final_free_dev); + STATUS("Overall log residual: initial = %e, final = %e\n", + init_log_dev, final_log_dev); + STATUS("Overall log free residual: initial = %e, final = %e\n", + init_free_log_dev, final_free_log_dev); check_rejection(crystals, n_crystals, full, max_B); normalise_scales(crystals, n_crystals); diff --git a/src/post-refinement.c b/src/post-refinement.c index d21108fd..1923ff4c 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -819,8 +819,8 @@ static void write_residual_graph(Crystal *cr, const RefList *full) } -static double do_scale_refine(Crystal *cr, const RefList *full, - PartialityModel pmodel, int verbose) +static void do_scale_refine(Crystal *cr, const RefList *full, + PartialityModel pmodel, int verbose) { int i, done; double old_dev; @@ -859,13 +859,11 @@ static double do_scale_refine(Crystal *cr, const RefList *full, STATUS("Final G=%.2f, B=%e\n", crystal_get_osf(cr), crystal_get_Bfac(cr)); } - - return old_dev; } -static double do_pr_refine(Crystal *cr, const RefList *full, - PartialityModel pmodel, int verbose) +static void do_pr_refine(Crystal *cr, const RefList *full, + PartialityModel pmodel, int verbose) { int i, done; double old_dev; @@ -917,8 +915,6 @@ static double do_pr_refine(Crystal *cr, const RefList *full, &csx, &csy, &csz); STATUS("Final asx = %e\n", asx); } - - return old_dev; } @@ -927,17 +923,17 @@ static struct prdata pr_refine(Crystal *cr, const RefList *full, { int verbose = 0; struct prdata prdata; - double final_dev; prdata.refined = 0; prdata.n_filtered = 0; prdata.initial_residual = residual(cr, full, 0, 0, NULL); prdata.initial_free_residual = residual(cr, full, 0, 1, NULL); - final_dev = prdata.initial_residual; + prdata.initial_log_residual = log_residual(cr, full, 0); + prdata.initial_free_log_residual = log_residual(cr, full, 1); if ( !no_scale ) { - final_dev = do_scale_refine(cr, full, pmodel, verbose); + do_scale_refine(cr, full, pmodel, verbose); } if ( verbose ) { @@ -945,16 +941,18 @@ static struct prdata pr_refine(Crystal *cr, const RefList *full, } if ( !no_pr ) { - final_dev = do_pr_refine(cr, full, pmodel, verbose); + do_pr_refine(cr, full, pmodel, verbose); } if ( crystal_get_user_flag(cr) == 0 ) { prdata.refined = 1; } + prdata.final_residual = residual(cr, full, 0, 0, NULL); prdata.final_free_residual = residual(cr, full, 0, 1, NULL); - prdata.final_residual = final_dev; + prdata.final_log_residual = log_residual(cr, full, 0); + prdata.final_free_log_residual = log_residual(cr, full, 1); - if ( isnan(final_dev) ) { + if ( isnan(prdata.final_residual) ) { STATUS("nan residual! Image serial %i\n", crystal_get_image(cr)->serial); } @@ -983,8 +981,12 @@ struct queue_args struct refine_args task_defaults; double initial_residual; double initial_free_residual; + double initial_log_residual; + double initial_free_log_residual; double final_residual; double final_free_residual; + double final_log_residual; + double final_free_log_residual; }; @@ -1016,16 +1018,32 @@ static void *get_image(void *vqargs) static void done_image(void *vqargs, void *task) { - struct queue_args *qargs = vqargs; + struct queue_args *qa = vqargs; struct refine_args *pargs = task; + struct prdata *prd = &pargs->prdata; + + qa->n_done++; + + if ( !isnan(prd->initial_residual) + && !isnan(prd->initial_free_residual) + && !isnan(prd->initial_log_residual) + && !isnan(prd->initial_free_log_residual) + && !isnan(prd->final_residual) + && !isnan(prd->final_free_residual) + && !isnan(prd->final_log_residual) + && !isnan(prd->final_free_log_residual) ) + { + qa->initial_residual += prd->initial_residual; + qa->initial_free_residual += prd->initial_free_residual; + qa->initial_log_residual += prd->initial_log_residual; + qa->initial_free_log_residual += prd->initial_free_log_residual; + qa->final_residual += prd->final_residual; + qa->final_free_residual += prd->final_free_residual; + qa->final_log_residual += prd->final_log_residual; + qa->final_free_log_residual += prd->final_free_log_residual; + } - qargs->n_done++; - qargs->initial_residual += pargs->prdata.initial_residual; - qargs->initial_free_residual += pargs->prdata.initial_free_residual; - qargs->final_residual += pargs->prdata.final_residual; - qargs->final_free_residual += pargs->prdata.final_free_residual; - - progress_bar(qargs->n_done, qargs->n_crystals, "Refining"); + progress_bar(qa->n_done, qa->n_crystals, "Refining"); free(task); } @@ -1034,7 +1052,9 @@ void refine_all(Crystal **crystals, int n_crystals, RefList *full, int nthreads, PartialityModel pmodel, int no_scale, int no_pr, double *initial_residual, double *initial_free_residual, - double *final_residual, double *final_free_residual) + double *initial_log_residual, double *initial_free_log_residual, + double *final_residual, double *final_free_residual, + double *final_log_residual, double *final_free_log_residual) { struct refine_args task_defaults; struct queue_args qargs; @@ -1054,8 +1074,12 @@ void refine_all(Crystal **crystals, int n_crystals, qargs.crystals = crystals; qargs.initial_residual = 0.0; qargs.initial_free_residual = 0.0; + qargs.initial_log_residual = 0.0; + qargs.initial_free_log_residual = 0.0; qargs.final_residual = 0.0; qargs.final_free_residual = 0.0; + qargs.final_log_residual = 0.0; + qargs.final_free_log_residual = 0.0; /* Don't have threads which are doing nothing */ if ( n_crystals < nthreads ) nthreads = n_crystals; @@ -1065,6 +1089,10 @@ void refine_all(Crystal **crystals, int n_crystals, *initial_residual = qargs.initial_residual; *initial_free_residual = qargs.initial_free_residual; + *initial_log_residual = qargs.initial_log_residual; + *initial_free_log_residual = qargs.initial_free_log_residual; *final_residual = qargs.final_residual; *final_free_residual = qargs.final_free_residual; + *final_log_residual = qargs.final_log_residual; + *final_free_log_residual = qargs.final_free_log_residual; } diff --git a/src/post-refinement.h b/src/post-refinement.h index fc72918f..f00a11ec 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -47,10 +47,18 @@ struct prdata { int refined; int n_filtered; + + /* Before refinement */ double initial_residual; double initial_free_residual; + double initial_log_residual; + double initial_free_log_residual; + + /* After refinement */ double final_residual; double final_free_residual; + double final_log_residual; + double final_free_log_residual; }; enum prflag @@ -69,8 +77,14 @@ extern const char *str_prflag(enum prflag flag); extern void refine_all(Crystal **crystals, int n_crystals, RefList *full, int nthreads, PartialityModel pmodel, int no_scale, int no_pr, - double *initial_residual, double *initial_free_residual, - double *final_residual, double *final_free_residual); + double *initial_residual, + double *initial_free_residual, + double *initial_log_residual, + double *initial_free_log_residual, + double *final_residual, + double *final_free_residual, + double *final_log_residual, + double *final_free_log_residual); /* Exported so it can be poked by tests/pr_p_gradient_check */ extern double gradient(Crystal *cr, int k, Reflection *refl, |