diff options
| -rw-r--r-- | src/partialator.c | 15 | ||||
| -rw-r--r-- | src/post-refinement.c | 23 | ||||
| -rw-r--r-- | src/post-refinement.h | 12 | ||||
| -rw-r--r-- | src/scaling-report.h | 3 |
4 files changed, 40 insertions, 13 deletions
diff --git a/src/partialator.c b/src/partialator.c index 6b29a74f..88912226 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -92,6 +92,7 @@ struct refine_args RefList *full; Crystal *crystal; PartialityModel pmodel; + struct prdata prdata; }; @@ -101,6 +102,7 @@ struct queue_args int n_done; Crystal **crystals; int n_crystals; + struct srdata *srdata; struct refine_args task_defaults; }; @@ -110,7 +112,7 @@ static void refine_image(void *task, int id) struct refine_args *pargs = task; Crystal *cr = pargs->crystal; - pr_refine(cr, pargs->full, pargs->pmodel); + pargs->prdata = pr_refine(cr, pargs->full, pargs->pmodel); } @@ -133,8 +135,10 @@ static void *get_image(void *vqargs) static void done_image(void *vqargs, void *task) { struct queue_args *qargs = vqargs; + struct refine_args *pargs = task; qargs->n_done++; + qargs->srdata->n_filtered += pargs->prdata.n_filtered; progress_bar(qargs->n_done, qargs->n_crystals, "Refining"); free(task); @@ -143,7 +147,8 @@ static void done_image(void *vqargs, void *task) static void refine_all(Crystal **crystals, int n_crystals, struct detector *det, - RefList *full, int nthreads, PartialityModel pmodel) + RefList *full, int nthreads, PartialityModel pmodel, + struct srdata *srdata) { struct refine_args task_defaults; struct queue_args qargs; @@ -161,6 +166,7 @@ static void refine_all(Crystal **crystals, int n_crystals, qargs.n_done = 0; qargs.n_crystals = n_crystals; qargs.crystals = crystals; + qargs.srdata = srdata; /* Don't have threads which are doing nothing */ if ( n_crystals < nthreads ) nthreads = n_crystals; @@ -629,11 +635,14 @@ int main(int argc, char *argv[]) STATUS("Post refinement cycle %i of %i\n", i+1, n_iter); + srdata.n_filtered = 0; + /* Refine the geometry of all patterns to get the best fit */ comp = (reference == NULL) ? full : reference; select_reflections_for_refinement(crystals, n_crystals, comp, have_reference); - refine_all(crystals, n_crystals, det, comp, nthreads, pmodel); + refine_all(crystals, n_crystals, det, comp, nthreads, pmodel, + &srdata); nobs = 0; for ( j=0; j<n_crystals; j++ ) { diff --git a/src/post-refinement.c b/src/post-refinement.c index c444ae31..628e33d2 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -314,7 +314,7 @@ static void apply_shift(Crystal *cr, int k, double shift) } -static void check_eigen(gsl_vector *e_val) +static int check_eigen(gsl_vector *e_val) { int i; double vmax, vmin; @@ -352,12 +352,16 @@ static void check_eigen(gsl_vector *e_val) if ( verbose ) { STATUS("Condition number: %e / %e = %5.2f\n", vmax, vmin, vmax/vmin); + + STATUS("%i out of %i eigenvalues filtered.\n", n_filt, n); } + + return n_filt; } -static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M) +static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *n_filt) { gsl_matrix *s_vec; gsl_vector *s_val; @@ -403,7 +407,7 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M) /* "SAS" is now "U" */ /* Filter the eigenvalues */ - check_eigen(s_val); + *n_filt = check_eigen(s_val); /* Calculate the vector SB, which is the RHS of the equation */ SB = gsl_vector_calloc(n); @@ -437,7 +441,7 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M) /* Perform one cycle of post refinement on 'image' against 'full' */ static double pr_iterate(Crystal *cr, const RefList *full, - PartialityModel pmodel) + PartialityModel pmodel, struct prdata *prdata) { gsl_matrix *M; gsl_vector *v; @@ -540,7 +544,7 @@ static double pr_iterate(Crystal *cr, const RefList *full, } max_shift = 0.0; - shifts = solve_svd(v, M); + shifts = solve_svd(v, M, &prdata->n_filtered); if ( shifts != NULL ) { for ( param=0; param<NUM_PARAMS; param++ ) { @@ -642,12 +646,14 @@ static void free_backup_crystal(Crystal *cr) } -void pr_refine(Crystal *cr, const RefList *full, PartialityModel pmodel) +struct prdata pr_refine(Crystal *cr, const RefList *full, + PartialityModel pmodel) { double max_shift, dev; int i; Crystal *backup; const int verbose = 0; + struct prdata prdata; if ( verbose ) { dev = guide_dev(cr, full); @@ -673,7 +679,8 @@ void pr_refine(Crystal *cr, const RefList *full, PartialityModel pmodel) cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); - max_shift = pr_iterate(cr, full, pmodel); + prdata.n_filtered = 0; + max_shift = pr_iterate(cr, full, pmodel, &prdata); update_partialities_2(cr, pmodel, &n_gained, &n_lost); @@ -694,4 +701,6 @@ void pr_refine(Crystal *cr, const RefList *full, PartialityModel pmodel) } while ( (max_shift > 50.0) && (i < MAX_CYCLES) ); free_backup_crystal(backup); + + return prdata; } diff --git a/src/post-refinement.h b/src/post-refinement.h index 7e13090b..b964975e 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -41,6 +41,7 @@ #include "utils.h" #include "crystal.h" #include "geometry.h" +#include "scaling-report.h" /* Refineable parameters. @@ -55,13 +56,20 @@ enum { REF_CSX, REF_CSY, REF_CSZ, - REF_DIV, NUM_PARAMS, + REF_DIV, REF_R, }; -extern void pr_refine(Crystal *cr, const RefList *full, PartialityModel pmodel); +struct prdata +{ + int n_filtered; +}; + + +extern struct prdata pr_refine(Crystal *cr, const RefList *full, + PartialityModel pmodel); /* Exported so it can be poked by tests/pr_gradient_check */ extern double p_gradient(Crystal *cr, int k, Reflection *refl, diff --git a/src/scaling-report.h b/src/scaling-report.h index 473e1759..1d430042 100644 --- a/src/scaling-report.h +++ b/src/scaling-report.h @@ -41,10 +41,11 @@ typedef struct _srcontext SRContext; /* Opaque */ /* Information is logged in this structure */ struct srdata { - int n_no_refine; /* Number that couldn't be refined */ Crystal **crystals; int n; RefList *full; + + int n_filtered; }; #ifdef HAVE_CAIRO |