diff options
Diffstat (limited to 'libcrystfel/src/utils.c')
| -rw-r--r-- | libcrystfel/src/utils.c | 145 |
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; |