More work on scaling

4 files changed, 273 insertions, 66 deletions

diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index 83c18e17..11306029 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -27,6 +27,10 @@
 #include "cell.h"
 
 
+/* Maximum number of iterations of NLSq scaling per macrocycle. */
+#define MAX_CYCLES (10)
+
+
 static void show_matrix_eqn(gsl_matrix *M, gsl_vector *v, int r)
 {
 	int i, j;
@@ -41,92 +45,276 @@ static void show_matrix_eqn(gsl_matrix *M, gsl_vector *v, int r)
 }
 
 
-/* Scale the stack of images */
-void scale_intensities(struct image *images, int n, const char *sym)
+static double gradient(int j, signed int h, signed int k, signed int l,
+                       struct image *images, int n, const char *sym)
+{
+	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++ ) {
+
+		image = &images[m];
+		spots = image->cpeaks;
+
+		for ( i=0; i<image->n_cpeaks; i++ ) {
+
+			signed int ha, ka, la;
+
+			if ( !spots[i].valid ) continue;
+			if ( spots[i].p < 0.1 ) continue;
+
+			get_asymm(spots[i].h, spots[i].k, spots[i].l,
+			          &ha, &ka, &la, sym);
+
+			if ( ha != h ) continue;
+			if ( ka != k ) continue;
+			if ( la != l ) continue;
+
+			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);
+
+		}
+
+	}
+
+	/* "This frame" terms */
+	image = &images[j];
+	spots = image->cpeaks;
+	for ( i=0; i<image->n_cpeaks; i++ ) {
+
+		signed int ha, ka, la;
+
+		if ( !spots[i].valid ) continue;
+		if ( spots[i].p < 0.1 ) continue;
+
+		get_asymm(spots[i].h, spots[i].k, spots[i].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);
+
+	}
+
+	return (num1*num1_this - num2*num2_this) / den;
+}
+
+
+static double iterate_scale(struct image *images, int n, const char *sym,
+                            double *I_full_old)
 {
-#if 0
 	gsl_matrix *M;
 	gsl_vector *v;
 	gsl_vector *shifts;
-	int j;
+	int m;
+	double max_shift;
+	int crossed = 0;
 
-	M = gsl_matrix_calloc(n, n);
-	v = gsl_vector_calloc(n);
+	M = gsl_matrix_calloc(n-1, n-1);
+	v = gsl_vector_calloc(n-1);
 
-	for ( j=0; j<n; j++ ) {
+	for ( m=0; m<n; m++ ) {  /* "Equation number": one equation per frame */
 
-		signed int hind, kind, lind;
-		signed int ha, ka, la;
-		double I_full, delta_I;
-		float I_partial;
-		float xc, yc;
+		int k;  /* Frame (scale factor) number */
 		int h;
-		struct image *image = &images[j];
+		int mcomp;
+		double vc_tot = 0.0;
+		struct image *image = &images[m];
 		struct cpeak *spots = image->cpeaks;
 
-		for ( h=0; h<image->n_cpeaks; h++ ) {
+		if ( m == crossed ) continue;
 
-			int g;
-			double v_c, gr, I_full;
+		/* Determine the "solution" vector component */
+		for ( h=0; h<image->n_cpeaks; h++ ) {
 
-			hind = spots[h].h;
-			kind = spots[h].k;
-			lind = spots[h].l;
+			double v_c;
+			double ic, ip;
+			signed int ha, ka, la;
 
-			/* Don't attempt to use spots with very small
-			 * partialities, since it won't be accurate. */
+			if ( !spots[h].valid ) continue;
 			if ( spots[h].p < 0.1 ) continue;
 
-			if ( integrate_peak(image, spots[h].x, spots[h].y,
-				            &xc, &yc, &I_partial, NULL, NULL, 1, 1) ) {
-				continue;
-			}
+			get_asymm(spots[h].h, spots[h].k, spots[h].l,
+			          &ha, &ka, &la, sym);
+			ic = lookup_intensity(I_full_old, ha, ka, la);
+
+			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;
+
+		}
+
+		mcomp = m;
+		if ( m > crossed ) mcomp--;
+		gsl_vector_set(v, mcomp, vc_tot);
 
-			get_asymm(hind, kind, lind, &ha, &ka, &la, sym);
-			I_full = lookup_intensity(i_full, ha, ka, la);
-			delta_I = I_partial - (spots[h].p * I_full / image->osf);
+		/* Now fill in the matrix components */
+		for ( k=0; k<n; k++ ) {
 
+			double val = 0.0;
+			int kcomp;
 
-			for ( g=0; g<NUM_PARAMS; g++ ) {
+			if ( k == crossed ) continue;
 
-				double M_curr, M_c;
+			for ( h=0; h<image->n_cpeaks; h++ ) {
 
-				M_curr = gsl_matrix_get(M, g, k);
+				signed int ha, ka, la;
+				double con;
+				double ic;
 
-				M_c = gradient(image, g, spots[h],
-				               image->profile_radius)
-				    * gradient(image, k, spots[h],
-				               image->profile_radius);
-				M_c *= pow(I_full, 2.0);
+				if ( !spots[h].valid ) continue;
+				if ( spots[h].p < 0.1 ) continue;
 
-				gsl_matrix_set(M, g, k, M_curr + M_c);
+				get_asymm(spots[h].h, spots[h].k, spots[h].l,
+				          &ha, &ka, &la, sym);
+				ic = lookup_intensity(I_full_old, ha, ka, la);
+
+				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;
 
 			}
 
-			gr = gradient(image, k, spots[h],
-			              image->profile_radius);
-			v_c = delta_I * I_full * gr;
-			gsl_vector_set(v, k, v_c);
+			kcomp = k;
+			if ( k > crossed ) kcomp--;
+			gsl_matrix_set(M, mcomp, kcomp, val);
 
 		}
 
+
 	}
-	show_matrix_eqn(M, v, NUM_PARAMS);
+	show_matrix_eqn(M, v, n-1);
 
-	shifts = gsl_vector_alloc(NUM_PARAMS);
+	shifts = gsl_vector_alloc(n-1);
 	gsl_linalg_HH_solve(M, v, shifts);
-	for ( param=0; param<NUM_PARAMS; param++ ) {
-		double shift = gsl_vector_get(shifts, param);
-		apply_shift(image, param, shift);
+	max_shift = 0.0;
+	for ( m=0; m<n-1; m++ ) {
+
+		double shift = gsl_vector_get(shifts, m);
+		int mimg;
+
+		mimg = m;
+		if ( mimg >= crossed ) mimg++;
+
+		images[mimg].osf += shift;
+
+		if ( shift > max_shift ) max_shift = shift;
+
 	}
+	gsl_vector_free(shifts);
 
 	gsl_matrix_free(M);
 	gsl_vector_free(v);
-	gsl_vector_free(shifts);
 
-	free(spots);
-	spots = find_intersections(image, image->indexed_cell, n, 0);
-	*pspots = spots;
-	return mean_partial_dev(image, spots, *n, sym, i_full, NULL);
-#endif
+	return max_shift;
+}
+
+
+static double *lsq_intensities(struct image *images, int n, const char *sym,
+                               ReflItemList *obs)
+{
+	double *I_full;
+	int i;
+
+	I_full = new_list_intensity();
+	for ( i=0; i<num_items(obs); i++ ) {
+
+		signed int h, k, l;
+		struct refl_item *it = get_item(obs, i);
+		double num = 0.0;
+		double den = 0.0;
+		int m;
+
+		get_asymm(it->h, it->k, it->l, &h, &k, &l, sym);
+
+		/* For each frame */
+		for ( m=0; m<n; m++ ) {
+
+			double G;
+			int a;
+
+			G = images[m].osf;
+
+			/* For each peak */
+			for ( a=0; a<images[m].n_cpeaks; a++ ) {
+
+				signed int ha, ka, la;
+
+				if ( !images[m].cpeaks[a].valid ) continue;
+				if ( images[m].cpeaks[a].p < 0.1 ) continue;
+
+				/* Correct reflection? */
+				get_asymm(images[m].cpeaks[a].h,
+				          images[m].cpeaks[a].k,
+				          images[m].cpeaks[a].l,
+				          &ha, &ka, &la, sym);
+				if ( ha != h ) continue;
+				if ( ka != k ) continue;
+				if ( la != l ) continue;
+
+				num += images[m].cpeaks[a].intensity
+				     * images[m].cpeaks[a].p * G;
+
+				den += pow(images[m].cpeaks[a].p, 2.0)
+				     * pow(G, 2.0);
+
+			}
+
+		}
+
+		set_intensity(I_full, h, k, l, num/den);
+
+	}
+
+	return I_full;
+}
+
+
+/* Scale the stack of images */
+double *scale_intensities(struct image *images, int n, const char *sym,
+                          ReflItemList *obs)
+{
+	int m;
+	double *I_full;
+	int i;
+	double max_shift;
+
+	/* 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);
+
+	/* Iterate */
+	i = 0;
+	do {
+
+		max_shift = iterate_scale(images, n, sym, I_full);
+		STATUS("Iteration %2i: max shift = %5.2f\n", i, max_shift);
+		free(I_full);
+		I_full = lsq_intensities(images, n, sym, obs);
+		i++;
+
+	} while ( (max_shift > 0.1) && (i < MAX_CYCLES) );
+
+	return I_full;
 }
diff --git a/src/hrs-scaling.h b/src/hrs-scaling.h
index 28343283..ae384aa8 100644
--- a/src/hrs-scaling.h
+++ b/src/hrs-scaling.h
@@ -20,7 +20,8 @@
 
 #include "image.h"
 
-extern void scale_intensities(struct image *image, int n, const char *sym);
+extern double *scale_intensities(struct image *image, int n, const char *sym,
+                                 ReflItemList *obs);
 
 
 #endif	/* HRS_SCALING_H */
diff --git a/src/image.h b/src/image.h
index 0d2f6c2e..078231a0 100644
--- a/src/image.h
+++ b/src/image.h
@@ -61,6 +61,7 @@ struct cpeak
 	signed int l;
 
 	double min_distance;
+	int valid;
 
 	/* Partiality */
 	double r1;  /* First excitation error */
@@ -72,6 +73,9 @@ struct cpeak
 	/* Location in image */
 	int x;
 	int y;
+
+	/* Intensity */
+	double intensity;
 };
 
 
diff --git a/src/partialator.c b/src/partialator.c
index c023fcde..b90c8c27 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -163,7 +163,7 @@ static void refine_all(struct image *images, int n_total_patterns,
 }
 
 
-static void integrate_image(struct image *image)
+static void integrate_image(struct image *image, ReflItemList *obs)
 {
 	struct cpeak *spots;
 	int j, n;
@@ -210,15 +210,22 @@ static void integrate_image(struct image *image)
 		if ( integrate_peak(image, spots[j].x, spots[j].y,
 		                    &xc, &yc, &i_partial, NULL, NULL,
 		                    1, 1, 0) ) {
+			spots[j].valid = 0;
 			continue;
 		}
 
+		spots[j].intensity = i_partial;
+		spots[j].valid = 1;
+
+		if ( !find_item(obs, h, k, l) ) add_item(obs, h, k, l);
+
 	}
+	image->cpeaks = spots;
+	image->n_cpeaks = n;
 
 	free(image->data);
 	if ( image->flags != NULL ) free(image->flags);
 	hdfile_close(hdfile);
-	image->cpeaks = spots;
 
 	/* Muppet proofing */
 	image->data = NULL;
@@ -239,7 +246,6 @@ int main(int argc, char *argv[])
 	int config_basename = 0;
 	int config_checkprefix = 1;
 	struct detector *det;
-	double *i_full;
 	unsigned int *cts;
 	ReflItemList *obs;
 	int i;
@@ -247,6 +253,7 @@ int main(int argc, char *argv[])
 	struct image *images;
 	int n_iter = 10;
 	struct beam_params *beam = NULL;
+	double *I_full;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -363,10 +370,6 @@ int main(int argc, char *argv[])
 		return 1;
 	}
 
-	/* Prepare for iteration */
-	i_full = new_list_intensity();
-	obs = new_items();
-
 	n_total_patterns = count_patterns(fh);
 	STATUS("There are %i patterns to process\n", n_total_patterns);
 
@@ -378,6 +381,7 @@ int main(int argc, char *argv[])
 
 	/* Fill in what we know about the images so far */
 	rewind(fh);
+	obs = new_items();
 	for ( i=0; i<n_total_patterns; i++ ) {
 
 		UnitCell *cell;
@@ -416,7 +420,7 @@ int main(int argc, char *argv[])
 		images[i].flags = NULL;
 
 		/* Get reflections from this image */
-		integrate_image(&images[i]);
+		integrate_image(&images[i], obs);
 
 		progress_bar(i, n_total_patterns-1, "Loading pattern data");
 
@@ -428,7 +432,7 @@ int main(int argc, char *argv[])
 
 	/* Make initial estimates */
 	STATUS("Performing initial scaling.\n");
-	scale_intensities(images, n_total_patterns, sym);
+	I_full = scale_intensities(images, n_total_patterns, sym, obs);
 
 	/* Iterate */
 	for ( i=0; i<n_iter; i++ ) {
@@ -454,21 +458,31 @@ int main(int argc, char *argv[])
 		}
 
 		/* Refine the geometry of all patterns to get the best fit */
-		refine_all(images, n_total_patterns, det, sym, obs, i_full,
+		refine_all(images, n_total_patterns, det, sym, obs, I_full,
 		           nthreads, fhg, fhp);
 
 		/* Re-estimate all the full intensities */
-		scale_intensities(images, n_total_patterns, sym);
+		free(I_full);
+		I_full = scale_intensities(images, n_total_patterns, sym, obs);
 
 		fclose(fhg);
 		fclose(fhp);
+
+	}
+
+	STATUS("Final scale factors:\n");
+	for ( i=0; i<n_total_patterns; i++ ) {
+		STATUS("%4i : %e\n", i, images[i].osf);
 	}
 
 	/* Output results */
-	write_reflections(outfile, obs, i_full, NULL, NULL, cts, NULL);
+	write_reflections(outfile, obs, I_full, NULL, NULL, cts, NULL);
 
 	/* Clean up */
-	free(i_full);
+	for ( i=0; i<n_total_patterns; i++ ) {
+		free(images[i].cpeaks);
+	}
+	free(I_full);
 	delete_items(obs);
 	free(sym);
 	free(outfile);