Move solve_svd() to utils

4 files changed, 149 insertions, 224 deletions

diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index e18f6681..a05f64f1 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -54,94 +54,6 @@
 #include "integration.h"
 
 
-static void check_eigen(gsl_vector *e_val)
-{
-	int i;
-	double vmax, vmin;
-	const int n = e_val->size;
-	const double max_condition = 1e6;
-	const int verbose = 0;
-
-	if ( verbose ) STATUS("Eigenvalues:\n");
-	vmin = +INFINITY;
-	vmax = 0.0;
-	for ( i=0; i<n; i++ ) {
-		double val = gsl_vector_get(e_val, i);
-		if ( verbose ) STATUS("%i: %e\n", i, val);
-		if ( val > vmax ) vmax = val;
-		if ( val < vmin ) vmin = val;
-	}
-
-	for ( i=0; i<n; i++ ) {
-		double val = gsl_vector_get(e_val, i);
-		if ( val < vmax/max_condition ) {
-			gsl_vector_set(e_val, i, 0.0);
-		}
-	}
-
-	vmin = +INFINITY;
-	vmax = 0.0;
-	for ( i=0; i<n; i++ ) {
-		double val = gsl_vector_get(e_val, i);
-		if ( val == 0.0 ) continue;
-		if ( val > vmax ) vmax = val;
-		if ( val < vmin ) vmin = val;
-	}
-	if ( verbose ) {
-		STATUS("Condition number: %e / %e = %5.2f\n",
-		       vmax, vmin, vmax/vmin);
-	}
-}
-
-
-static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *Mp)
-{
-	gsl_matrix *s_vec;
-	gsl_vector *s_val;
-	int err, n;
-	gsl_vector *shifts;
-	gsl_matrix *M;
-
-	n = v->size;
-	if ( v->size != Mp->size1 ) return NULL;
-	if ( v->size != Mp->size2 ) return NULL;
-
-	M = gsl_matrix_alloc(n, n);
-	if ( M == NULL ) return NULL;
-	gsl_matrix_memcpy(M, Mp);
-
-	s_val = gsl_vector_calloc(n);
-	s_vec = gsl_matrix_calloc(n, n);
-
-	err = gsl_linalg_SV_decomp_jacobi(M, s_vec, s_val);
-	if ( err ) {
-		ERROR("SVD failed: %s\n", gsl_strerror(err));
-		gsl_matrix_free(s_vec);
-		gsl_vector_free(s_val);
-		return NULL;
-	}
-	/* "M" is now "U" */
-
-	check_eigen(s_val);
-
-	shifts = gsl_vector_calloc(n);
-	err = gsl_linalg_SV_solve(M, s_vec, s_val, v, shifts);
-	if ( err ) {
-		ERROR("Matrix solution failed: %s\n", gsl_strerror(err));
-		gsl_matrix_free(s_vec);
-		gsl_vector_free(s_val);
-		gsl_vector_free(shifts);
-		return NULL;
-	}
-
-	gsl_matrix_free(s_vec);
-	gsl_vector_free(s_val);
-	gsl_matrix_free(M);
-
-	return shifts;
-}
-
-
 enum boxmask_val
 {
 	BM_IG,  /* "Soft" ignore */
@@ -445,8 +357,8 @@ static void fit_bg(struct intcontext *ic, struct peak_box *bx)
 	}
 	}
 
-	/* SVD is massive overkill here */
-	ans = solve_svd(v, bx->bgm);
+	/* FIXME: SVD is massive overkill here */
+	ans = solve_svd(v, bx->bgm, NULL, 0);
 	gsl_vector_free(v);
 
 	bx->a = gsl_vector_get(ans, 0);
diff --git a/libcrystfel/src/utils.c b/libcrystfel/src/utils.c
index f8f93046..643e0158 100644
--- a/libcrystfel/src/utils.c
+++ b/libcrystfel/src/utils.c
@@ -40,6 +40,8 @@
 #include <gsl/gsl_matrix.h>
 #include <gsl/gsl_vector.h>
 #include <gsl/gsl_blas.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_eigen.h>
 
 #include "utils.h"
 #include "image.h"
@@ -108,6 +110,149 @@ void show_matrix(gsl_matrix *M)
 }
 
 
+static int check_eigen(gsl_vector *e_val, int verbose)
+{
+	int i;
+	double vmax, vmin;
+	const int n = e_val->size;
+	const double max_condition = 1e6;
+	int n_filt = 0;
+
+	if ( verbose ) STATUS("Eigenvalues:\n");
+	vmin = +INFINITY;
+	vmax = 0.0;
+	for ( i=0; i<n; i++ ) {
+		double val = gsl_vector_get(e_val, i);
+		if ( verbose ) STATUS("%i: %e\n", i, val);
+		if ( val > vmax ) vmax = val;
+		if ( val < vmin ) vmin = val;
+	}
+
+	for ( i=0; i<n; i++ ) {
+		double val = gsl_vector_get(e_val, i);
+		if ( val < vmax/max_condition ) {
+			gsl_vector_set(e_val, i, 0.0);
+			n_filt++;
+		}
+	}
+
+	vmin = +INFINITY;
+	vmax = 0.0;
+	for ( i=0; i<n; i++ ) {
+		double val = gsl_vector_get(e_val, i);
+		if ( val == 0.0 ) continue;
+		if ( val > vmax ) vmax = val;
+		if ( val < vmin ) vmin = 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;
+}
+
+
+/**
+ * solve_svd:
+ * v: a gsl_vector
+ * M: a gsl_matrix
+ * n_filt: pointer to store the number of filtered eigenvalues
+ * verbose: flag for verbosity on the terminal
+ *
+ * Solves the matrix equation M.x = v, returning x.
+ * Performs rescaling and eigenvalue filtering.
+ **/
+gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *n_filt, int verbose)
+{
+	gsl_matrix *s_vec;
+	gsl_vector *s_val;
+	int err, n;
+	gsl_vector *shifts;
+	gsl_vector *SB;
+	gsl_vector *SinvX;
+	gsl_matrix *S;  /* rescaling matrix due to Bricogne */
+	gsl_matrix *AS;
+	gsl_matrix *SAS;
+	int i;
+	gsl_matrix *SAS_copy;
+
+	n = v->size;
+	if ( v->size != M->size1 ) return NULL;
+	if ( v->size != M->size2 ) return NULL;
+
+	/* Calculate the rescaling matrix S */
+	S = gsl_matrix_calloc(n, n);
+	for ( i=0; i<n; i++ ) {
+		double sii = pow(gsl_matrix_get(M, i, i), -0.5);
+		gsl_matrix_set(S, i, i, sii);
+	}
+
+	/* Calculate the matrix SAS, which we will be (not) inverting */
+	AS = gsl_matrix_calloc(n, n);
+	SAS = gsl_matrix_calloc(n, n);
+	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, M, S, 0.0, AS);
+	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);
+	err = gsl_linalg_SV_decomp_jacobi(SAS, s_vec, s_val);
+	if ( err ) {
+		if ( verbose ) ERROR("SVD failed: %s\n", gsl_strerror(err));
+		gsl_matrix_free(s_vec);
+		gsl_vector_free(s_val);
+		gsl_matrix_free(SAS);
+		gsl_matrix_free(S);
+		return NULL;
+	}
+	/* "SAS" is now "U" */
+
+	/* Filter the eigenvalues */
+	*n_filt = check_eigen(s_val, verbose);
+
+	gsl_matrix_free(SAS_copy);
+
+	/* Solve the equation SAS.SinvX = SB */
+	SinvX = gsl_vector_calloc(n);
+	err = gsl_linalg_SV_solve(SAS, s_vec, s_val, SB, SinvX);
+	gsl_vector_free(SB);
+	gsl_matrix_free(SAS);
+	gsl_matrix_free(s_vec);
+	gsl_vector_free(s_val);
+
+	if ( err ) {
+		ERROR("Matrix solution failed: %s\n", gsl_strerror(err));
+		gsl_matrix_free(S);
+		gsl_vector_free(SinvX);
+		return NULL;
+	}
+
+	/* Calculate S.SinvX to get X, the shifts */
+	shifts = gsl_vector_calloc(n);
+	gsl_blas_dgemv(CblasNoTrans, 1.0, S, SinvX, 0.0, shifts);
+
+	gsl_matrix_free(S);
+	gsl_vector_free(SinvX);
+
+	return shifts;
+}
+
+
 size_t notrail(char *s)
 {
 	ssize_t i;
diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h
index 2f2009d6..4955f875 100644
--- a/libcrystfel/src/utils.h
+++ b/libcrystfel/src/utils.h
@@ -103,6 +103,8 @@ extern struct rvec quat_rot(struct rvec q, struct quaternion z);
 
 extern void show_matrix_eqn(gsl_matrix *M, gsl_vector *v);
 extern void show_matrix(gsl_matrix *M);
+extern gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M, int *n_filt,
+                            int verbose);
 extern size_t notrail(char *s);
 extern void chomp(char *s);
 
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 4dae5096..96c0de12 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -362,140 +362,6 @@ static void apply_shift(Crystal *cr, int k, double shift)
 }
 
 
-static int check_eigen(gsl_vector *e_val, int verbose)
-{
-	int i;
-	double vmax, vmin;
-	const int n = e_val->size;
-	const double max_condition = 1e6;
-	int n_filt = 0;
-
-	if ( verbose ) STATUS("Eigenvalues:\n");
-	vmin = +INFINITY;
-	vmax = 0.0;
-	for ( i=0; i<n; i++ ) {
-		double val = gsl_vector_get(e_val, i);
-		if ( verbose ) STATUS("%i: %e\n", i, val);
-		if ( val > vmax ) vmax = val;
-		if ( val < vmin ) vmin = val;
-	}
-
-	for ( i=0; i<n; i++ ) {
-		double val = gsl_vector_get(e_val, i);
-		if ( val < vmax/max_condition ) {
-			gsl_vector_set(e_val, i, 0.0);
-			n_filt++;
-		}
-	}
-
-	vmin = +INFINITY;
-	vmax = 0.0;
-	for ( i=0; i<n; i++ ) {
-		double val = gsl_vector_get(e_val, i);
-		if ( val == 0.0 ) continue;
-		if ( val > vmax ) vmax = val;
-		if ( val < vmin ) vmin = 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, int *n_filt,
-                             int verbose)
-{
-	gsl_matrix *s_vec;
-	gsl_vector *s_val;
-	int err, n;
-	gsl_vector *shifts;
-	gsl_vector *SB;
-	gsl_vector *SinvX;
-	gsl_matrix *S;  /* rescaling matrix due to Bricogne */
-	gsl_matrix *AS;
-	gsl_matrix *SAS;
-	int i;
-	gsl_matrix *SAS_copy;
-
-	n = v->size;
-	if ( v->size != M->size1 ) return NULL;
-	if ( v->size != M->size2 ) return NULL;
-
-	/* Calculate the rescaling matrix S */
-	S = gsl_matrix_calloc(n, n);
-	for ( i=0; i<n; i++ ) {
-		double sii = pow(gsl_matrix_get(M, i, i), -0.5);
-		gsl_matrix_set(S, i, i, sii);
-	}
-
-	/* Calculate the matrix SAS, which we will be (not) inverting */
-	AS = gsl_matrix_calloc(n, n);
-	SAS = gsl_matrix_calloc(n, n);
-	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, M, S, 0.0, AS);
-	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);
-	err = gsl_linalg_SV_decomp_jacobi(SAS, s_vec, s_val);
-	if ( err ) {
-		if ( verbose ) ERROR("SVD failed: %s\n", gsl_strerror(err));
-		gsl_matrix_free(s_vec);
-		gsl_vector_free(s_val);
-		gsl_matrix_free(SAS);
-		gsl_matrix_free(S);
-		return NULL;
-	}
-	/* "SAS" is now "U" */
-
-	/* Filter the eigenvalues */
-	*n_filt = check_eigen(s_val, verbose);
-
-	gsl_matrix_free(SAS_copy);
-
-	/* Solve the equation SAS.SinvX = SB */
-	SinvX = gsl_vector_calloc(n);
-	err = gsl_linalg_SV_solve(SAS, s_vec, s_val, SB, SinvX);
-	gsl_vector_free(SB);
-	gsl_matrix_free(SAS);
-	gsl_matrix_free(s_vec);
-	gsl_vector_free(s_val);
-
-	if ( err ) {
-		ERROR("Matrix solution failed: %s\n", gsl_strerror(err));
-		gsl_matrix_free(S);
-		gsl_vector_free(SinvX);
-		return NULL;
-	}
-
-	/* Calculate S.SinvX to get X, the shifts */
-	shifts = gsl_vector_calloc(n);
-	gsl_blas_dgemv(CblasNoTrans, 1.0, S, SinvX, 0.0, shifts);
-
-	gsl_matrix_free(S);
-	gsl_vector_free(SinvX);
-
-	return shifts;
-}
-
-
 /* Perform one cycle of post refinement on 'image' against 'full' */
 static double pr_iterate(Crystal *cr, const RefList *full,
                          PartialityModel pmodel, int *n_filtered)