diff options
Diffstat (limited to 'src/hrs-scaling.c')
-rw-r--r-- | src/hrs-scaling.c | 136 |
1 files changed, 89 insertions, 47 deletions
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index bde7ad98..9b5005d8 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -135,8 +135,62 @@ static double s_vh(struct image *images, int n, } +static gsl_vector *solve_by_eigenvalue_filtration(gsl_vector *v, gsl_matrix *M) +{ + gsl_eigen_symmv_workspace *work; + gsl_vector *e_val; + gsl_matrix *e_vec; + int val, n, frame; + gsl_vector *shifts; + + n = v->size; + if ( v->size != M->size1 ) return NULL; + if ( v->size != M->size2 ) return NULL; + + /* Diagonalise */ + work = gsl_eigen_symmv_alloc(n); + e_val = gsl_vector_alloc(n); + e_vec = gsl_matrix_alloc(n, n); + val = gsl_eigen_symmv(M, e_val, e_vec, work); + if ( val ) { + ERROR("Couldn't diagonalise matrix.\n"); + return NULL; + } + gsl_eigen_symmv_free(work); + val = gsl_eigen_symmv_sort(e_val, e_vec, GSL_EIGEN_SORT_ABS_DESC); + + /* Rotate the "solution vector" */ + gsl_vector *rprime; + rprime = gsl_vector_alloc(n); + val = gsl_blas_dgemv(CblasTrans, 1.0, e_vec, v, 0.0, rprime); + + /* Solve (now easy) */ + gsl_vector *sprime; + sprime = gsl_vector_alloc(n); + for ( frame=0; frame<n-1; frame++ ) { + double num, den; + num = gsl_vector_get(rprime, frame); + den = gsl_vector_get(e_val, frame); + gsl_vector_set(sprime, frame, num/den); + } + gsl_vector_set(sprime, n-1, 0.0); /* Set last shift to zero */ + + /* Rotate back */ + shifts = gsl_vector_alloc(n); + val = gsl_blas_dgemv(CblasNoTrans, 1.0, e_vec, sprime, 0.0, shifts); + + gsl_vector_free(sprime); + gsl_vector_free(rprime); + gsl_matrix_free(e_vec); + gsl_vector_free(e_val); + + return shifts; +} + + static double iterate_scale(struct image *images, int n, - ReflItemList *obs, const char *sym, char *cref) + ReflItemList *obs, const char *sym, char *cref, + double *reference) { gsl_matrix *M; gsl_vector *v; @@ -197,16 +251,26 @@ static double iterate_scale(struct image *images, int n, int b; /* Frame (scale factor) number */ struct image *image_a = &images[a]; - double vc, Ih, uh, vh, rha, vha, uha; + double vc, Ih, uh, rha, vha, uha; double vval; /* Determine the "solution" vector component */ uha = uha_arr[a]; vha = vha_arr[a]; uh = lookup_intensity(uh_arr, h, k, l); - vh = lookup_intensity(vh_arr, h, k, l); - Ih = vh / uh; - if ( isnan(Ih) ) Ih = 0.0; /* 0 / 0 = 0, not NaN */ + + + if ( !reference ) { + double vh; + vh = lookup_intensity(vh_arr, h, k, l); + Ih = vh / uh; + /* 0 / 0 = 0, not NaN */ + if ( isnan(Ih) ) Ih = 0.0; + } else { + /* Look up by asymmetric indices */ + Ih = lookup_intensity(reference, h, k, l); + } + rha = vha - image_a->osf * uha * Ih; vc = Ih * rha; @@ -220,6 +284,9 @@ static double iterate_scale(struct image *images, int n, /* Matrix is symmetric */ if ( b > a ) continue; + /* Off-diagonals zero if reference available */ + if ( (reference != NULL) && (a!=b) ) continue; + /* Zero if no common reflections */ if ( cref[a+n*b] != 0 ) continue; @@ -252,43 +319,24 @@ static double iterate_scale(struct image *images, int n, free(uha_arr); free(vha_arr); - /* Fox and Holmes method */ - gsl_eigen_symmv_workspace *work; - gsl_vector *e_val; - gsl_matrix *e_vec; - int val; - - /* Diagonalise */ - work = gsl_eigen_symmv_alloc(n); - e_val = gsl_vector_alloc(n); - e_vec = gsl_matrix_alloc(n, n); - val = gsl_eigen_symmv(M, e_val, e_vec, work); - if ( val ) { - ERROR("Couldn't diagonalise matrix.\n"); - return 0.0; + show_matrix_eqn(M, v, n); + + if ( reference == NULL ) { + shifts = solve_by_eigenvalue_filtration(v, M); + } else { + shifts = gsl_vector_alloc(n); + for ( frame=0; frame<n-1; frame++ ) { + double num, den; + num = gsl_vector_get(v, frame); + den = gsl_matrix_get(M, frame, frame); + gsl_vector_set(shifts, frame, num/den); + } } - gsl_eigen_symmv_free(work); - val = gsl_eigen_symmv_sort(e_val, e_vec, GSL_EIGEN_SORT_ABS_DESC); - /* Rotate the "solution vector" */ - gsl_vector *rprime; - rprime = gsl_vector_alloc(n); - val = gsl_blas_dgemv(CblasTrans, 1.0, e_vec, v, 0.0, rprime); - - /* Solve (now easy) */ - gsl_vector *sprime; - sprime = gsl_vector_alloc(n); - for ( frame=0; frame<n-1; frame++ ) { - double num, den; - num = gsl_vector_get(rprime, frame); - den = gsl_vector_get(e_val, frame); - gsl_vector_set(sprime, frame, num/den); - } - gsl_vector_set(sprime, n-1, 0.0); /* Set last shift to zero */ + gsl_matrix_free(M); + gsl_vector_free(v); - /* Rotate back */ - shifts = gsl_vector_alloc(n); - val = gsl_blas_dgemv(CblasNoTrans, 1.0, e_vec, sprime, 0.0, shifts); + if ( shifts == NULL ) return 0.0; /* Apply shifts */ max_shift = 0.0; @@ -305,12 +353,6 @@ static double iterate_scale(struct image *images, int n, } gsl_vector_free(shifts); - gsl_vector_free(sprime); - gsl_vector_free(rprime); - gsl_matrix_free(e_vec); - gsl_vector_free(e_val); - gsl_matrix_free(M); - gsl_vector_free(v); return max_shift; } @@ -394,7 +436,7 @@ static void normalise_osfs(struct image *images, int n) /* Scale the stack of images */ RefList *scale_intensities(struct image *images, int n, const char *sym, - ReflItemList *obs, char *cref) + ReflItemList *obs, char *cref, double *reference) { int m; RefList *full; @@ -408,7 +450,7 @@ RefList *scale_intensities(struct image *images, int n, const char *sym, i = 0; do { - max_shift = iterate_scale(images, n, obs, sym, cref); + max_shift = iterate_scale(images, n, obs, sym, cref, reference); STATUS("Scaling iteration %2i: max shift = %5.2f\n", i, max_shift); i++; |