From d05eecd314661fad45db4fb48df179f861d77d24 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Mon, 12 Aug 2013 10:41:30 +0200
Subject: partialator: Debugging stuff

---
 src/partialator.c     | 25 +++++++++++++++++++++++++
 src/post-refinement.c | 22 ++++++++++++++++------
 src/post-refinement.h |  2 +-
 src/scaling-report.c  |  3 +++
 4 files changed, 45 insertions(+), 7 deletions(-)

diff --git a/src/partialator.c b/src/partialator.c
index 88912226..ad6d9a49 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -309,6 +309,16 @@ static void display_progress(int n_images, int n_crystals)
 }
 
 
+static const char *str_flags(Crystal *cr)
+{
+	if ( crystal_get_user_flag(cr) ) {
+		return "N";
+	}
+
+	return "-";
+}
+
+
 int main(int argc, char *argv[])
 {
 	int c;
@@ -671,6 +681,21 @@ int main(int argc, char *argv[])
 	/* Output results */
 	write_reflist(outfile, full);
 
+	/* Dump parameters */
+	FILE *fh;
+	fh = fopen("partialator.params", "w");
+	if ( fh == NULL ) {
+		ERROR("Couldn't open partialator.params!\n");
+	} else {
+		for ( i=0; i<n_crystals; i++ ) {
+			fprintf(fh, "%4i %5.2f %8.5e %s\n", i,
+			        crystal_get_osf(crystals[i]),
+			        crystal_get_image(crystals[i])->div,
+			        str_flags(crystals[i]));
+		}
+		fclose(fh);
+	}
+
 	/* Clean up */
 	for ( i=0; i<n_crystals; i++ ) {
 		reflist_free(crystal_get_reflections(crystals[i]));
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 94c8a816..88c609c1 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -314,13 +314,12 @@ static void apply_shift(Crystal *cr, int k, double shift)
 }
 
 
-static int check_eigen(gsl_vector *e_val)
+static int check_eigen(gsl_matrix *M, gsl_vector *e_val, int verbose)
 {
 	int i;
 	double vmax, vmin;
 	const int n = e_val->size;
 	const double max_condition = 1e6;
-	const int verbose = 0;
 	int n_filt = 0;
 
 	if ( verbose ) STATUS("Eigenvalues:\n");
@@ -374,6 +373,7 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *n_filt,
 	gsl_matrix *AS;
 	gsl_matrix *SAS;
 	int i;
+	gsl_matrix *SAS_copy;
 
 	n = v->size;
 	if ( v->size != M->size1 ) return NULL;
@@ -393,6 +393,18 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *n_filt,
 	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, S, AS, 0.0, SAS);
 	gsl_matrix_free(AS);
 
+	/* Calculate the vector SB, which is the RHS of the equation */
+	SB = gsl_vector_calloc(n);
+	gsl_blas_dgemv(CblasNoTrans, 1.0, S, v, 0.0, SB);
+
+	if ( verbose ) {
+		STATUS("The equation after rescaling:\n");
+		show_matrix_eqn(SAS, SB);
+	}
+
+	SAS_copy = gsl_matrix_alloc(SAS->size1, SAS->size2);
+	gsl_matrix_memcpy(SAS_copy, SAS);
+
 	/* Do the SVD */
 	s_val = gsl_vector_calloc(n);
 	s_vec = gsl_matrix_calloc(n, n);
@@ -408,11 +420,9 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *n_filt,
 	/* "SAS" is now "U" */
 
 	/* Filter the eigenvalues */
-	*n_filt = check_eigen(s_val);
+	*n_filt = check_eigen(SAS_copy, s_val, verbose);
 
-	/* Calculate the vector SB, which is the RHS of the equation */
-	SB = gsl_vector_calloc(n);
-	gsl_blas_dgemv(CblasNoTrans, 1.0, S, v, 0.0, SB);
+	gsl_matrix_free(SAS_copy);
 
 	/* Solve the equation SAS.SinvX = SB */
 	SinvX = gsl_vector_calloc(n);
diff --git a/src/post-refinement.h b/src/post-refinement.h
index b964975e..63b6ff66 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -56,8 +56,8 @@ enum {
 	REF_CSX,
 	REF_CSY,
 	REF_CSZ,
-	NUM_PARAMS,
 	REF_DIV,
+	NUM_PARAMS,
 	REF_R,
 };
 
diff --git a/src/scaling-report.c b/src/scaling-report.c
index acb007ae..17268489 100644
--- a/src/scaling-report.c
+++ b/src/scaling-report.c
@@ -889,6 +889,9 @@ void sr_iteration(SRContext *sr, int iteration, struct srdata *d)
 		cairo_restore(sr->cr);
 
 	}
+
+	STATUS("%i filtered total, %5.2f filtered per crystal\n",
+	       d->n_filtered, (double)d->n_filtered / d->n);
 }
 
 
-- 
cgit v1.2.3