aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src/utils.c
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel/src/utils.c')
-rw-r--r--libcrystfel/src/utils.c145
1 files changed, 145 insertions, 0 deletions
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;