aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2013-08-01 17:33:23 +0200
committerThomas White <taw@physics.org>2013-08-01 17:33:23 +0200
commit4c0971b30a4982c8b4a3be4b6599639ecb1b6695 (patch)
tree48b419284a4273b87f96c97c1cdee70f8dee4cab
parentd001bcca749215e41a79a0de32e4cc049ace8b86 (diff)
Count filtered eigenvalues
-rw-r--r--src/partialator.c15
-rw-r--r--src/post-refinement.c23
-rw-r--r--src/post-refinement.h12
-rw-r--r--src/scaling-report.h3
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