From 0431728cc5efd3321b1c60f97a830dd525cf04c8 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Thu, 25 Jul 2013 12:04:01 +0200
Subject: partialator: Simplify scaling

---
 src/hrs-scaling.c | 38 +++++++++++++++++++++-----------------
 src/partialator.c |  3 ++-
 2 files changed, 23 insertions(+), 18 deletions(-)

(limited to 'src')

diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index c249ef46..47c4301e 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -48,10 +48,7 @@
 
 
 /* Maximum number of iterations of scaling per macrocycle. */
-#define MAX_CYCLES (50)
-
-/* ESD of restraint driving scale factors to unity */
-#define SCALING_RESTRAINT (1.0)
+#define MAX_CYCLES (10)
 
 
 struct scale_queue_args
@@ -110,9 +107,8 @@ static void run_scale_job(void *vwargs, int cookie)
 	      refl = next_refl(refl, iter) )
 	{
 		signed int h, k, l;
-		double Ih, Ihl, esd;
+		double Ih, Ihl, esd, corr;
 		Reflection *r;
-		double corr;
 
 		if ( !get_scalable(refl) ) continue;
 
@@ -151,19 +147,18 @@ static void run_scale_job(void *vwargs, int cookie)
 		Ihl = get_intensity(refl) / corr;
 		esd = get_esd_intensity(refl) / corr;
 
-		num += Ih * (Ihl/G) / pow(esd/G, 2.0);
-		den += pow(Ih, 2.0)/pow(esd/G, 2.0);
+		num += Ih * Ihl;
+		den += Ihl * Ihl;
 
 	}
 
-	//num += image->osf / pow(SCALING_RESTRAINT, 2.0);
-	//den += pow(image->osf, 2.0)/pow(SCALING_RESTRAINT, 2.0);
-
 	g = num / den;
 	if ( !isnan(g) && !isinf(g) ) {
-		crystal_set_osf(cr, g*G);
+		crystal_set_osf(cr, g);
+		wargs->shift = fabs((g/G)-1.0);
+	} else {
+		wargs->shift = 0.0;
 	}
-	wargs->shift = fabs(g-1.0);
 }
 
 
@@ -298,8 +293,8 @@ static void run_merge_job(void *vwargs, int cookie)
 		Ihl = get_intensity(refl) / corr;
 		esd = get_esd_intensity(refl) / corr;
 
-		num += (Ihl/G) / pow(esd/G, 2.0);
-		den += 1.0 / pow(esd/G, 2.0);
+		num += Ihl * G;
+		den += 1.0;
 		red++;
 
 		set_temp1(f, num);
@@ -518,6 +513,7 @@ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref,
 	i = 0;
 	do {
 
+		int j;
 		RefList *reference;
 
 		/* Refine against reference or current "full" estimates */
@@ -529,8 +525,12 @@ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref,
 
 		max_corr = iterate_scale(crystals, n, reference, n_threads,
 		                         pmodel);
-		//STATUS("Scaling iteration %2i: max correction = %5.2f\n",
-		//       i+1, max_corr);
+		STATUS("Scaling iteration %2i: max correction = %5.2f\n",
+		       i+1, max_corr);
+		for ( j=0; j<10; j++ ) {
+			printf(" %5.2f", crystal_get_osf(crystals[j]));
+		}
+		printf("\n");
 
 		/* No reference -> generate list for next iteration */
 		if ( gref == NULL ) {
@@ -542,6 +542,10 @@ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref,
 
 	} while ( (max_corr > 0.01) && (i < MAX_CYCLES) );
 
+	if ( i == MAX_CYCLES ) {
+		ERROR("Warning: Scaling did not converge.\n");
+	}
+
 	if ( gref != NULL ) {
 		full = lsq_intensities(crystals, n, n_threads, pmodel);
 	} /* else we already did it */
diff --git a/src/partialator.c b/src/partialator.c
index 73046e0a..07d7437c 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -190,7 +190,7 @@ static int select_scalable_reflections(RefList *list, RefList *reference)
 		if ( get_partiality(refl) < 0.1 ) sc = 0;
 		v = fabs(get_intensity(refl));
 		esd = get_esd_intensity(refl);
-		if ( v < 0.5*esd ) sc = 0;
+		//if ( v < 0.5*esd ) sc = 0;
 
 		/* If we are scaling against a reference set, we additionally
 		 * require that this reflection is in the reference list. */
@@ -497,6 +497,7 @@ int main(int argc, char *argv[])
 	n_crystals = 0;
 	images = NULL;
 	crystals = NULL;
+
 	do {
 
 		RefList *as;
-- 
cgit v1.2.3