diff options
| -rw-r--r-- | src/partialator.c | 38 | ||||
| -rw-r--r-- | src/post-refinement.c | 112 | ||||
| -rw-r--r-- | src/post-refinement.h | 32 |
3 files changed, 64 insertions, 118 deletions
diff --git a/src/partialator.c b/src/partialator.c index de08df34..81c5ab37 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -585,6 +585,20 @@ static void normalise_scales(Crystal **crystals, int n_crystals) } +static void show_all_residuals(Crystal **crystals, int n_crystals, + RefList *full) +{ + double dev, free_dev, log_dev, free_log_dev; + + all_residuals(crystals, n_crystals, full, + &dev, &free_dev, &log_dev, &free_log_dev); + STATUS("Residuals:" + " linear linear free log log free\n"); + STATUS(" "); + STATUS("%15e %15e %15e %15e\n", dev, free_dev, log_dev, free_log_dev); +} + + int main(int argc, char *argv[]) { int c; @@ -990,29 +1004,15 @@ int main(int argc, char *argv[]) /* Iterate */ 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; + show_all_residuals(crystals, n_crystals, full); 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, max_B, - &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); + no_scale, no_pr, max_B); + + show_all_residuals(crystals, n_crystals, full); check_rejection(crystals, n_crystals, full); normalise_scales(crystals, n_crystals); @@ -1026,6 +1026,8 @@ int main(int argc, char *argv[]) } + show_all_residuals(crystals, n_crystals, full); + /* Output results */ STATUS("Writing overall results to %s\n", outfile); write_reflist_2(outfile, full, sym); diff --git a/src/post-refinement.c b/src/post-refinement.c index 9d1aa7a2..818911df 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -54,6 +54,12 @@ /* Maximum number of iterations of NLSq to do for each image per macrocycle. */ #define MAX_CYCLES (10) +struct prdata +{ + int refined; + int n_filtered; +}; + const char *str_prflag(enum prflag flag) { switch ( flag ) { @@ -774,6 +780,38 @@ static double residual(Crystal *cr, const RefList *full, int verbose, 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, 0, NULL); + free_r = residual(crystals[i], full, 0, 1, NULL); + log_r = log_residual(crystals[i], full, 0); + free_log_r = log_residual(crystals[i], full, 1); + + 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; @@ -928,11 +966,6 @@ static struct prdata pr_refine(Crystal *cr, const RefList *full, 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); - prdata.initial_log_residual = log_residual(cr, full, 0); - prdata.initial_free_log_residual = log_residual(cr, full, 1); - if ( !no_scale ) { do_scale_refine(cr, full, pmodel, verbose); } @@ -954,15 +987,6 @@ static struct prdata pr_refine(Crystal *cr, const RefList *full, 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_log_residual = log_residual(cr, full, 0); - prdata.final_free_log_residual = log_residual(cr, full, 1); - - if ( isnan(prdata.final_residual) ) { - STATUS("nan residual! Image serial %i\n", - crystal_get_image(cr)->serial); - } return prdata; } @@ -987,14 +1011,6 @@ struct queue_args Crystal **crystals; int n_crystals; 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; }; @@ -1027,30 +1043,9 @@ static void *get_image(void *vqargs) static void done_image(void *vqargs, void *task) { 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; - } - progress_bar(qa->n_done, qa->n_crystals, "Refining"); free(task); } @@ -1058,11 +1053,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, - 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) + int no_scale, int no_pr, double max_B) { struct refine_args task_defaults; struct queue_args qargs; @@ -1075,41 +1066,16 @@ void refine_all(Crystal **crystals, int n_crystals, task_defaults.no_scale = no_scale; task_defaults.no_pr = no_pr; task_defaults.max_B = max_B; - task_defaults.prdata.initial_residual = 0.0; - task_defaults.prdata.initial_free_residual = 0.0; - task_defaults.prdata.initial_log_residual = 0.0; - task_defaults.prdata.initial_free_log_residual = 0.0; - task_defaults.prdata.final_residual = 0.0; - task_defaults.prdata.final_free_residual = 0.0; - task_defaults.prdata.final_log_residual = 0.0; - task_defaults.prdata.final_free_log_residual = 0.0; qargs.task_defaults = task_defaults; qargs.n_started = 0; qargs.n_done = 0; qargs.n_crystals = 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; run_threads(nthreads, refine_image, get_image, done_image, &qargs, n_crystals, 0, 0, 0); - - *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 9c0f9984..b3509501 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -43,24 +43,6 @@ #include "geometry.h" -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 { PRFLAG_OK = 0, @@ -76,18 +58,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 max_B, - 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); + int no_scale, int no_pr, double max_B); /* Exported so it can be poked by tests/pr_p_gradient_check */ extern double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel); +extern void all_residuals(Crystal **crystals, int n_crystals, RefList *full, + double *residual, double *free_residual, + double *log_residual, double *free_log_residual); + #endif /* POST_REFINEMENT_H */ |