aboutsummaryrefslogtreecommitdiff
path: root/src/hrs-scaling.c
diff options
context:
space:
mode:
authorThomas White <taw@bitwiz.org.uk>2010-12-17 15:56:01 -0800
committerThomas White <taw@physics.org>2012-02-22 15:27:10 +0100
commitcdaa57c6db66fe243ed00730a18395a2a01e8712 (patch)
tree2ed9de46978cb61464a9f91dfcba7b1ea7ee7254 /src/hrs-scaling.c
parent4dccca06469f97b7b13abab8bc9cd1dac71a04c1 (diff)
More work on scaling
Diffstat (limited to 'src/hrs-scaling.c')
-rw-r--r--src/hrs-scaling.c254
1 files changed, 148 insertions, 106 deletions
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index 11306029..80894131 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -28,7 +28,7 @@
/* Maximum number of iterations of NLSq scaling per macrocycle. */
-#define MAX_CYCLES (10)
+#define MAX_CYCLES (30)
static void show_matrix_eqn(gsl_matrix *M, gsl_vector *v, int r)
@@ -45,156 +45,167 @@ static void show_matrix_eqn(gsl_matrix *M, gsl_vector *v, int r)
}
-static double gradient(int j, signed int h, signed int k, signed int l,
- struct image *images, int n, const char *sym)
+static void sum_GI(struct image *images, int n, const char *sym,
+ signed int hat, signed int kat, signed int lat,
+ double *sigma_GI, double *sigma_Gsq)
{
- struct image *image;
- struct cpeak *spots;
- double num1 = 0.0;
- double num2 = 0.0;
- double den = 0.0;
- double num1_this = 0.0;
- double num2_this = 0.0;
- int m, i;
-
- /* "Everything" terms */
- for ( m=0; m<n; m++ ) {
+ int k;
- image = &images[m];
- spots = image->cpeaks;
+ *sigma_GI = 0.0;
+ *sigma_Gsq = 0.0;
+ for ( k=0; k<n; k++ ) {
- for ( i=0; i<image->n_cpeaks; i++ ) {
+ int hi;
+ struct image *image = &images[k];
+ struct cpeak *spots = images->cpeaks;
+ int found = 0;
- signed int ha, ka, la;
+ for ( hi=0; hi<image->n_cpeaks; hi++ ) {
- if ( !spots[i].valid ) continue;
- if ( spots[i].p < 0.1 ) continue;
+ double ic;
+ signed int ha, ka, la;
- get_asymm(spots[i].h, spots[i].k, spots[i].l,
- &ha, &ka, &la, sym);
+ if ( !spots[hi].valid ) continue;
+ if ( spots[hi].p < 0.1 ) continue;
+ get_asymm(spots[hi].h, spots[hi].k, spots[hi].l,
+ &ha, &ka, &la, sym);
+ if ( ha != hat ) continue;
+ if ( ka != kat ) continue;
+ if ( la != lat ) continue;
- if ( ha != h ) continue;
- if ( ka != k ) continue;
- if ( la != l ) continue;
+ ic = spots[hi].intensity / spots[hi].p;
- num1 += pow(image->osf, 2.0) * pow(spots[i].p, 2.0);
- num2 += image->osf * spots[i].intensity * spots[i].p;
- den += pow(image->osf, 2.0) * pow(spots[i].p, 2.0);
+ *sigma_GI += ic * image->osf;
+ found = 1;
}
+ if ( found ) {
+ *sigma_Gsq += pow(image->osf, 2.0);
+ }
}
+}
- /* "This frame" terms */
- image = &images[j];
- spots = image->cpeaks;
- for ( i=0; i<image->n_cpeaks; i++ ) {
- signed int ha, ka, la;
+static double find_occurrances(struct image *image, const char *sym,
+ signed int h, signed int k, signed int l)
+{
+ double Ihl = 0.0;
+ int find;
+ struct cpeak *spots = image->cpeaks;
- if ( !spots[i].valid ) continue;
- if ( spots[i].p < 0.1 ) continue;
+ for ( find=0; find<image->n_cpeaks; find++ ) {
- get_asymm(spots[i].h, spots[i].k, spots[i].l,
- &ha, &ka, &la, sym);
+ signed int ha, ka, la;
+ if ( !spots[find].valid ) continue;
+ if ( spots[find].p < 0.1 ) continue;
+ get_asymm(spots[find].h, spots[find].k,
+ spots[find].l, &ha, &ka, &la, sym);
if ( ha != h ) continue;
if ( ka != k ) continue;
if ( la != l ) continue;
- num1_this += spots[i].intensity * spots[i].p;
- num2_this += pow(spots[i].p, 2.0);
+ Ihl += spots[find].intensity / spots[find].p;
}
- return (num1*num1_this - num2*num2_this) / den;
+ return Ihl;
}
-static double iterate_scale(struct image *images, int n, const char *sym,
- double *I_full_old)
+static double iterate_scale(struct image *images, int n,
+ ReflItemList *obs, const char *sym)
{
gsl_matrix *M;
gsl_vector *v;
gsl_vector *shifts;
- int m;
+ int l;
double max_shift;
int crossed = 0;
+ int n_ref;
M = gsl_matrix_calloc(n-1, n-1);
v = gsl_vector_calloc(n-1);
+ n_ref = num_items(obs);
- for ( m=0; m<n; m++ ) { /* "Equation number": one equation per frame */
+ for ( l=0; l<n; l++ ) { /* "Equation number": one equation per frame */
- int k; /* Frame (scale factor) number */
+ int m; /* Frame (scale factor) number */
int h;
- int mcomp;
+ int lcomp;
double vc_tot = 0.0;
- struct image *image = &images[m];
- struct cpeak *spots = image->cpeaks;
+ struct image *imagel = &images[l];
- if ( m == crossed ) continue;
+ if ( l == crossed ) continue;
/* Determine the "solution" vector component */
- for ( h=0; h<image->n_cpeaks; h++ ) {
+ for ( h=0; h<n_ref; h++ ) {
- double v_c;
- double ic, ip;
- signed int ha, ka, la;
+ double sigma_GI, sigma_Gsq;
+ double vc;
+ double Ihl;
+ struct refl_item *it = get_item(obs, h);
- if ( !spots[h].valid ) continue;
- if ( spots[h].p < 0.1 ) continue;
+ sum_GI(images, n, sym, it->h, it->k, it->l,
+ &sigma_GI, &sigma_Gsq);
- get_asymm(spots[h].h, spots[h].k, spots[h].l,
- &ha, &ka, &la, sym);
- ic = lookup_intensity(I_full_old, ha, ka, la);
+ /* Add up symmetric equivalents within the pattern */
+ Ihl = find_occurrances(imagel, sym,
+ it->h, it->k, it->l);
- v_c = ip - (spots[h].p * image->osf * ic);
- v_c *= pow(spots[h].p, 2.0);
- v_c *= pow(image->osf, 2.0);
- v_c *= ic;
- v_c *= gradient(m, ha, ka, la, images, n, sym);
- vc_tot += v_c;
+ vc = Ihl * sigma_GI / sigma_Gsq;
+ vc -= imagel->osf * pow(sigma_GI, 2.0) / sigma_Gsq;
+
+ vc_tot += vc;
}
- mcomp = m;
- if ( m > crossed ) mcomp--;
- gsl_vector_set(v, mcomp, vc_tot);
+ lcomp = l;
+ if ( l > crossed ) lcomp--;
+ gsl_vector_set(v, lcomp, vc_tot);
/* Now fill in the matrix components */
- for ( k=0; k<n; k++ ) {
+ for ( m=0; m<n; m++ ) {
- double val = 0.0;
- int kcomp;
+ double mc_tot = 0.0;
+ int mcomp;
+ struct image *imagem = &images[m];
- if ( k == crossed ) continue;
+ if ( m == crossed ) continue;
- for ( h=0; h<image->n_cpeaks; h++ ) {
+ for ( h=0; h<n_ref; h++ ) {
- signed int ha, ka, la;
- double con;
- double ic;
+ double mc = 0.0;
+ double Ihl, Ihm;
+ struct refl_item *it = get_item(obs, h);
+ double sigma_GI, sigma_Gsq;
- if ( !spots[h].valid ) continue;
- if ( spots[h].p < 0.1 ) continue;
+ sum_GI(images, n, sym, it->h, it->k, it->l,
+ &sigma_GI, &sigma_Gsq);
- get_asymm(spots[h].h, spots[h].k, spots[h].l,
- &ha, &ka, &la, sym);
- ic = lookup_intensity(I_full_old, ha, ka, la);
+ if ( l == m ) {
+ mc += pow(sigma_GI, 2.0)
+ / pow(sigma_Gsq, 2.0);
+ }
+
+ Ihl = find_occurrances(imagel, sym,
+ it->h, it->k, it->l);
+ Ihm = find_occurrances(imagem, sym,
+ it->h, it->k, it->l);
- con = -pow(spots[h].p, 3.0);
- con *= pow(image->osf, 3.0);
- con *= ic;
- con *= gradient(m, ha, ka, la, images, n, sym);
- con *= gradient(k, ha, ka, la, images, n, sym);
- val += con;
+ mc += Ihl * Ihm / sigma_Gsq;
+
+ mc += (sigma_GI / pow(sigma_Gsq, 2.0) )
+ * ( imagel->osf*Ihm + imagem->osf * Ihl);
+
+ mc_tot += mc;
}
- kcomp = k;
- if ( k > crossed ) kcomp--;
- gsl_matrix_set(M, mcomp, kcomp, val);
+ mcomp = m;
+ if ( m > crossed ) mcomp--;
+ gsl_matrix_set(M, lcomp, mcomp, mc_tot);
}
@@ -205,17 +216,19 @@ static double iterate_scale(struct image *images, int n, const char *sym,
shifts = gsl_vector_alloc(n-1);
gsl_linalg_HH_solve(M, v, shifts);
max_shift = 0.0;
- for ( m=0; m<n-1; m++ ) {
+ for ( l=0; l<n-1; l++ ) {
- double shift = gsl_vector_get(shifts, m);
- int mimg;
+ double shift = gsl_vector_get(shifts, l);
+ int limg;
- mimg = m;
- if ( mimg >= crossed ) mimg++;
+ limg = l;
+ if ( limg >= crossed ) limg++;
- images[mimg].osf += shift;
+ images[limg].osf += shift;
- if ( shift > max_shift ) max_shift = shift;
+ if ( fabs(shift) > fabs(max_shift) ) {
+ max_shift = fabs(shift);
+ }
}
gsl_vector_free(shifts);
@@ -227,8 +240,8 @@ static double iterate_scale(struct image *images, int n, const char *sym,
}
-static double *lsq_intensities(struct image *images, int n, const char *sym,
- ReflItemList *obs)
+static double *lsq_intensities(struct image *images, int n,
+ ReflItemList *obs, const char *sym)
{
double *I_full;
int i;
@@ -287,6 +300,23 @@ static double *lsq_intensities(struct image *images, int n, const char *sym,
}
+static void normalise_osfs(struct image *images, int n)
+{
+ int m;
+ double tot = 0.0;
+ double mean;
+
+ for ( m=0; m<n; m++ ) {
+ tot += images[m].osf;
+ }
+ mean = tot / (double)n;
+
+ for ( m=0; m<n; m++ ) {
+ images[m].osf /= mean;
+ }
+}
+
+
/* Scale the stack of images */
double *scale_intensities(struct image *images, int n, const char *sym,
ReflItemList *obs)
@@ -298,23 +328,35 @@ double *scale_intensities(struct image *images, int n, const char *sym,
/* Start with all scale factors equal */
for ( m=0; m<n; m++ ) {
- images[m].osf = 1.0;
- }
- /* Calculate LSQ intensities using these scale factors */
- I_full = lsq_intensities(images, n, sym, obs);
+ double tot = 0.0;
+ int j;
+
+ for ( j=0; j<images[m].n_cpeaks; j++ ) {
+ tot += images[m].cpeaks[j].intensity
+ / images[m].cpeaks[j].p;
+ }
+
+ images[m].osf = tot;
+
+ }
+ normalise_osfs(images, n);
/* Iterate */
i = 0;
do {
- max_shift = iterate_scale(images, n, sym, I_full);
+ max_shift = iterate_scale(images, n, obs, sym);
STATUS("Iteration %2i: max shift = %5.2f\n", i, max_shift);
- free(I_full);
- I_full = lsq_intensities(images, n, sym, obs);
i++;
+ normalise_osfs(images, n);
+
+ } while ( (max_shift > 0.01) && (i < MAX_CYCLES) );
- } while ( (max_shift > 0.1) && (i < MAX_CYCLES) );
+ for ( m=0; m<n; m++ ) {
+ images[m].osf /= images[0].osf;
+ }
+ I_full = lsq_intensities(images, n, obs, sym);
return I_full;
}