Count filtered eigenvalues

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