Work in progress on post refinement

3 files changed, 55 insertions, 28 deletions

diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index cb0f8cde..fcf7efae 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -401,7 +401,8 @@ RefList *scale_intensities(struct image *images, int n, const char *sym,
 	do {
 
 		max_shift = iterate_scale(images, n, obs, sym, cref);
-		STATUS("Iteration %2i: max shift = %5.2f\n", i, max_shift);
+		STATUS("Scaling iteration %2i: max shift = %5.2f\n",
+		       i, max_shift);
 		i++;
 		normalise_osfs(images, n);
 
diff --git a/src/partial_sim.c b/src/partial_sim.c
index a9abb55b..2c2493e9 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -39,7 +39,7 @@ static void mess_up_cell(UnitCell *cell)
 	double cx, cy, cz;
 
 	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
-	ax += 0.005*ax;
+	ax += 0.008*ax;
 	cell_set_cartesian(cell, ax, ay, az, bx, by, bz, cx, cy, cz);
 }
 
@@ -255,7 +255,7 @@ int main(int argc, char *argv[])
 
 	/* Alter the cell by a tiny amount */
 	image.filename = "(simulated 2)";
-	new = rotate_cell(cell, deg2rad(0.001), deg2rad(0.0), 0.0);
+	new = rotate_cell(cell, deg2rad(1.0), deg2rad(0.0), 0.0);
 	cell_free(image.indexed_cell);
 	image.indexed_cell = new;
 
@@ -264,7 +264,7 @@ int main(int argc, char *argv[])
 	calculate_partials(image.reflections, 0.5, full, sym);
 
 	/* Give a slightly incorrect cell in the stream */
-	//mess_up_cell(image.indexed_cell);
+	mess_up_cell(image.indexed_cell);
 	write_chunk(ofh, &image, STREAM_INTEGRATED);
 	reflist_free(image.reflections);
 
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 7d851c17..20b8eeee 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -271,57 +271,63 @@ static double pr_iterate(struct image *image, const RefList *full,
 	      refl != NULL;
 	      refl = next_refl(refl, iter) ) {
 
-		signed int hind, kind, lind;
 		signed int ha, ka, la;
 		double I_full, delta_I;
 		double I_partial;
 		int k;
 		double p;
 		Reflection *match;
-
-		get_indices(refl, &hind, &kind, &lind);
+		double gradients[NUM_PARAMS];
 
 		if ( !get_scalable(refl) ) continue;
 
-		/* Actual measurement of this reflection from this pattern? */
-		I_partial = get_intensity(refl);
-
-		get_asymm(hind, kind, lind, &ha, &ka, &la, sym);
+		/* Find the full version */
+		get_indices(refl, &ha, &ka, &la);
 		match = find_refl(full, ha, ka, la);
-		assert(match != NULL);
+		assert(match != NULL);  /* Never happens because all scalable
+		                         * reflections had their LSQ intensities
+		                         * calculated in lsq_intensities(). */
 		I_full = get_intensity(match);
+
+		/* Actual measurement of this reflection from this pattern? */
+		I_partial = get_intensity(refl);
 		p = get_partiality(refl);
-		delta_I = I_partial - (p * I_full / image->osf);
+		delta_I = I_partial - (p * image->osf * I_full);
+
+		/* Calculate all gradients for this reflection */
+		for ( k=0; k<NUM_PARAMS; k++ ) {
+			double gr;
+			gr = gradient(image, k, refl, image->profile_radius);
+			gradients[k] = gr;
+		}
 
 		for ( k=0; k<NUM_PARAMS; k++ ) {
 
 			int g;
-			double v_c, gr;
+			double v_c, v_curr;
+			double gr;
 
 			for ( g=0; g<NUM_PARAMS; g++ ) {
 
-				double M_curr, M_c;
+				double M_c, M_curr;
 
-				M_curr = gsl_matrix_get(M, g, k);
-
-				M_c = gradient(image, g, refl,
-				               image->profile_radius)
-				    * gradient(image, k, refl,
-				               image->profile_radius);
+				M_c = gradients[g] * gradients[k];
 				M_c *= pow(I_full, 2.0);
 
+				M_curr = gsl_matrix_get(M, g, k);
 				gsl_matrix_set(M, g, k, M_curr + M_c);
 
 			}
 
 			gr = gradient(image, k, refl, image->profile_radius);
 			v_c = delta_I * I_full * gr;
-			gsl_vector_set(v, k, v_c);
+			v_curr = gsl_vector_get(v, k);
+			gsl_vector_set(v, k, v_curr + v_c);
 
 		}
 
 	}
-	//show_matrix_eqn(M, v, NUM_PARAMS);
+	show_matrix_eqn(M, v, NUM_PARAMS);
 
 	shifts = gsl_vector_alloc(NUM_PARAMS);
 	gsl_linalg_HH_solve(M, v, shifts);
@@ -381,15 +387,15 @@ static double mean_partial_dev(struct image *image,
 }
 
 
-void pr_refine(struct image *image, const RefList *full, const char *sym)
+static void plot_curve(struct image *image, const RefList *full,
+                       const char *sym)
 {
-	double max_shift;
-	int i;
 	double ax, ay, az;
 	double bx, by, bz;
 	double cx, cy, cz;
 	UnitCell *cell = image->indexed_cell;
 	double shval, origval;
+	int i;
 
 	cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
 	shval = 0.001*ax;
@@ -411,12 +417,32 @@ void pr_refine(struct image *image, const RefList *full, const char *sym)
 		STATUS("%i %e %e\n", i, ax, dev);
 
 	}
-	return;
+}
+
+
+void pr_refine(struct image *image, const RefList *full, const char *sym)
+{
+	double max_shift;
+	int i;
+
+	/* FIXME: This is for debugging */
+	//plot_curve(image, full, sym);
+	//return;
 
 	i = 0;
 	do {
+
+		double dev;
+
 		max_shift = pr_iterate(image, full, sym);
-		STATUS("Iteration %2i: max shift = %5.2f\n", i, max_shift);
+		update_partialities_and_asymm(image, sym,
+		                              NULL, NULL, NULL, NULL);
+
+		dev = mean_partial_dev(image, full, sym);
+		STATUS("PR Iteration %2i: max shift = %5.2f dev = %5.2f\n",
+		       i, max_shift, dev);
+
 		i++;
+
 	} while ( (max_shift > 0.01) && (i < MAX_CYCLES) );
 }