From 9225115a425ced041133e849f3816178967e3307 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 13 Mar 2015 11:06:28 +0100 Subject: Move solve_svd() to utils --- libcrystfel/src/integration.c | 92 +-------------------------- libcrystfel/src/utils.c | 145 ++++++++++++++++++++++++++++++++++++++++++ libcrystfel/src/utils.h | 2 + 3 files changed, 149 insertions(+), 90 deletions(-) (limited to 'libcrystfel/src') 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 vmax ) vmax = val; - if ( val < vmin ) vmin = val; - } - - for ( i=0; i 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 #include #include +#include +#include #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 vmax ) vmax = val; + if ( val < vmin ) vmin = val; + } + + for ( i=0; i 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; isize1, 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); -- cgit v1.2.3