aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/hrs-scaling.c136
-rw-r--r--src/hrs-scaling.h3
-rw-r--r--src/partialator.c19
3 files changed, 104 insertions, 54 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++;
diff --git a/src/hrs-scaling.h b/src/hrs-scaling.h
index d65129b8..f65b7988 100644
--- a/src/hrs-scaling.h
+++ b/src/hrs-scaling.h
@@ -21,7 +21,8 @@
#include "image.h"
extern RefList *scale_intensities(struct image *images, int n, const char *sym,
- ReflItemList *obs, char *cref);
+ ReflItemList *obs, char *cref,
+ double *reference);
extern char *find_common_reflections(struct image *images, int n);
diff --git a/src/partialator.c b/src/partialator.c
index 64546a7b..e3eccafb 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -171,7 +171,7 @@ int main(int argc, char *argv[])
char *cref;
int n_usable_patterns = 0;
char *reference_file = NULL;
- RefList *reference;
+ double *reference = NULL;
/* Long options */
const struct option longopts[] = {
@@ -230,7 +230,7 @@ int main(int argc, char *argv[])
break;
case 1 :
- reference = strdup(optarg);
+ reference_file = strdup(optarg);
break;
case 0 :
@@ -278,9 +278,15 @@ int main(int argc, char *argv[])
}
if ( reference_file != NULL ) {
- reference = read_reflections(reference_file);
+ RefList *list;
+ RefList *symmed;
+ list = read_reflections(reference_file);
free(reference_file);
- if ( reference == NULL ) return 1;
+ if ( list == NULL ) return 1;
+ symmed = asymmetric_indices(list, sym);
+ reflist_free(list);
+ reference = intensities_from_list(symmed);
+ reflist_free(symmed);
}
n_total_patterns = count_patterns(fh);
@@ -360,7 +366,7 @@ int main(int argc, char *argv[])
/* Make initial estimates */
STATUS("Performing initial scaling.\n");
full = scale_intensities(images, n_usable_patterns, sym,
- scalable, cref);
+ scalable, cref, reference);
for ( i=0; i<num_items(scalable); i++ ) {
Reflection *f;
@@ -425,7 +431,7 @@ int main(int argc, char *argv[])
/* Re-estimate all the full intensities */
reflist_free(full);
full = scale_intensities(images, n_usable_patterns,
- sym, scalable, cref);
+ sym, scalable, cref, reference);
fclose(fhg);
fclose(fhp);
@@ -451,6 +457,7 @@ int main(int argc, char *argv[])
free_detector_geometry(det);
free(beam);
free(cref);
+ free(reference);
for ( i=0; i<n_usable_patterns; i++ ) {
cell_free(images[i].indexed_cell);
free(images[i].filename);
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-02 09:12:58 +0000