Merge branch 'tom/logmerge'

6 files changed, 35 insertions, 17 deletions

diff --git a/src/merge.c b/src/merge.c
index 9d36a359..1dec13fd 100644
--- a/src/merge.c
+++ b/src/merge.c
@@ -63,6 +63,7 @@ struct merge_queue_args
 	double push_res;
 	int use_weak;
 	long long int n_reflections;
+	int ln_merge;
 };
 
 
@@ -168,6 +169,7 @@ static void run_merge_job(void *vwargs, int cookie)
 	Crystal *cr = wargs->crystal;
 	RefList *full = wargs->qargs->full;
 	double push_res = wargs->qargs->push_res;
+	int ln_merge = wargs->qargs->ln_merge;
 	Reflection *refl;
 	RefListIterator *iter;
 	double G, B;
@@ -193,7 +195,7 @@ static void run_merge_job(void *vwargs, int cookie)
 
 		if ( get_partiality(refl) < MIN_PART_MERGE ) continue;
 
-		if ( !wargs->qargs->use_weak ) {
+		if ( !wargs->qargs->use_weak || ln_merge ) {
 
 			if (get_intensity(refl) < 3.0*get_esd_intensity(refl)) {
 				continue;
@@ -235,7 +237,8 @@ static void run_merge_job(void *vwargs, int cookie)
 
 		/* Running mean and variance calculation */
 		temp = w + sumweight;
-		delta = get_intensity(refl)*corr - mean;
+		if (ln_merge) delta = log(get_intensity(refl)*corr) - mean;
+		else delta = get_intensity(refl)*corr - mean;
 		R = delta * w / temp;
 		set_intensity(f, mean + R);
 		set_temp2(f, M2 + sumweight * delta * R);
@@ -272,7 +275,7 @@ static void finalise_merge_job(void *vqargs, void *vwargs)
 
 RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
                            int min_meas,
-                           double push_res, int use_weak)
+                           double push_res, int use_weak, int ln_merge)
 {
 	RefList *full;
 	RefList *full2;
@@ -290,6 +293,7 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
 	qargs.push_res = push_res;
 	qargs.use_weak = use_weak;
 	qargs.n_reflections = 0;
+	qargs.ln_merge = ln_merge;
 	pthread_rwlock_init(&qargs.full_lock, NULL);
 
 	run_threads(n_threads, run_merge_job, create_merge_job,
@@ -305,10 +309,23 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
 	      refl != NULL;
 	      refl = next_refl(refl, iter) )
 	{
+
+		/* Correct for averaging log of intensities*/
+		if (ln_merge){
+			double ln_I, ln_temp2;
+
+			ln_temp2 = get_temp2(refl);
+			set_temp2(refl, exp(ln_temp2));
+
+			ln_I = get_intensity(refl);
+			set_intensity(refl, exp(ln_I));
+		}
+
 		double var;
 		int red;
 
 		red = get_redundancy(refl);
+		//TODO is this still correct for log averaging?
 		var = get_temp2(refl) / get_temp1(refl);
 		set_esd_intensity(refl, sqrt(var)/sqrt(red));
 
diff --git a/src/merge.h b/src/merge.h
index b33afdea..9bc54126 100644
--- a/src/merge.h
+++ b/src/merge.h
@@ -40,7 +40,8 @@
 #include "geometry.h"
 
 extern RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
-                                  int min_meas, double push_res, int use_weak);
+                                  int min_meas, double push_res, int use_weak,
+                                  int ln_merge);
 
 extern double residual(Crystal *cr, const RefList *full, int free,
                        int *pn_used, const char *filename);
diff --git a/src/partialator.c b/src/partialator.c
index 78b7173b..86431e4c 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -190,7 +190,7 @@ static void write_split(Crystal **crystals, int n_crystals, const char *outfile,
 	}
 	snprintf(tmp, 1024, "%s1", outfile);
 	split = merge_intensities(crystals1, n_crystals1, nthreads,
-		                  min_measurements, push_res, 1);
+		                  min_measurements, push_res, 1, 0);
 
 	if ( split == NULL ) {
 		ERROR("Not enough crystals for two way split!\n");
@@ -203,7 +203,7 @@ static void write_split(Crystal **crystals, int n_crystals, const char *outfile,
 	reflist_free(split);
 	snprintf(tmp, 1024, "%s2", outfile);
 	split = merge_intensities(crystals2, n_crystals2, nthreads,
-		                  min_measurements, push_res, 1);
+		                  min_measurements, push_res, 1, 0);
 	STATUS("and %s\n", tmp);
 	write_reflist_2(tmp, split, sym);
 	free_contribs(split);
@@ -291,7 +291,7 @@ static void write_custom_split(struct custom_split *csplit, int dsn,
 	STATUS("Writing dataset '%s' to %s (%i crystals)\n",
 	       csplit->dataset_names[dsn], tmp, n_crystalsn);
 	split = merge_intensities(crystalsn, n_crystalsn, nthreads,
-		                  min_measurements, push_res, 1);
+		                  min_measurements, push_res, 1, 0);
 	write_reflist_2(tmp, split, sym);
 	free_contribs(split);
 	reflist_free(split);
@@ -1412,7 +1412,7 @@ int main(int argc, char *argv[])
 			scale_all(crystals, n_crystals, nthreads, scaleflags);
 		}
 		full = merge_intensities(crystals, n_crystals, nthreads,
-		                         min_measurements, push_res, 1);
+		                         min_measurements, push_res, 1, 0);
 	} else {
 		full = reference;
 	}
@@ -1447,7 +1447,7 @@ int main(int argc, char *argv[])
 			}
 			full = merge_intensities(crystals, n_crystals, nthreads,
 			                         min_measurements,
-			                         push_res, 1);
+			                         push_res, 1, 0);
 		} /* else full still equals reference */
 
 		check_rejection(crystals, n_crystals, full, max_B, no_deltacchalf);
@@ -1495,10 +1495,10 @@ int main(int argc, char *argv[])
 		}
 		full = merge_intensities(crystals, n_crystals, nthreads,
 		                         min_measurements,
-		                         push_res, 1);
+		                         push_res, 1, 0);
 	} else {
 		full = merge_intensities(crystals, n_crystals, nthreads,
-		                         min_measurements, push_res, 1);
+		                         min_measurements, push_res, 1, 0);
 	}
 
 	/* Write final figures of merit (no rejection any more) */
diff --git a/src/scaling.c b/src/scaling.c
index 68456900..35481982 100644
--- a/src/scaling.c
+++ b/src/scaling.c
@@ -145,7 +145,7 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads, int scaleflags)
 		double bef_res;
 
 		full = merge_intensities(crystals, n_crystals, nthreads,
-		                         2, INFINITY, 0);
+		                         2, INFINITY, 0, 1);
 		old_res = new_res;
 		bef_res = total_log_r(crystals, n_crystals, full, NULL);
 
diff --git a/tests/partialator_merge_check_2 b/tests/partialator_merge_check_2
index a8ff3ddd..9938f64f 100755
--- a/tests/partialator_merge_check_2
+++ b/tests/partialator_merge_check_2
@@ -44,8 +44,8 @@ cat > partialator_merge_check_2_ans.hkl << EOF
 CrystFEL reflection list version 2.0
 Symmetry: 1
    h    k    l          I    phase   sigma(I)   nmeas
-   1    0    0     150.00        -       0.00       2
-  19    0    0      75.00        -       0.00       2
+   1    0    0     141.42        -       0.00       2
+  19    0    0      70.71        -       0.00       2
 EOF
 
 $PARTIALATOR -i partialator_merge_check_2.stream \
diff --git a/tests/partialator_merge_check_3 b/tests/partialator_merge_check_3
index 933422b0..2d237b02 100755
--- a/tests/partialator_merge_check_3
+++ b/tests/partialator_merge_check_3
@@ -45,9 +45,9 @@ cat > partialator_merge_check_3_ans.hkl << EOF
 CrystFEL reflection list version 2.0
 Symmetry: 4
    h    k    l          I    phase   sigma(I)   nmeas
-   1    0    0      75.00        -       0.00       2
-   2    0    0     150.00        -       0.00       2
-  20    0    0      75.00        -       0.00       2
+   1    0    0      70.71        -       0.00       2
+   2    0    0     141.42        -       0.00       2
+  20    0    0      70.71        -       0.00       2
 EOF
 
 $PARTIALATOR -i partialator_merge_check_3.stream \