1 files changed, 336 insertions, 124 deletions

diff --git a/src/partialator.c b/src/partialator.c
index 2540521d..e152db7c 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -3,11 +3,11 @@
  *
  * Scaling and post refinement for coherent nanocrystallography
  *
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2010-2015 Thomas White <taw@physics.org>
+ *   2010-2018 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -40,6 +40,7 @@
 #include <assert.h>
 #include <pthread.h>
 #include <gsl/gsl_errno.h>
+#include <sys/stat.h>
 
 #include <image.h>
 #include <utils.h>
@@ -189,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,
-		                  pmodel, min_measurements, push_res, 1);
+		                  min_measurements, push_res, 1);
 
 	if ( split == NULL ) {
 		ERROR("Not enough crystals for two way split!\n");
@@ -201,7 +202,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,
-		                  pmodel, min_measurements, push_res, 1);
+		                  min_measurements, push_res, 1);
 	STATUS("and %s\n", tmp);
 	write_reflist_2(tmp, split, sym);
 	reflist_free(split);
@@ -288,7 +289,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,
-		                  pmodel, min_measurements, push_res, 1);
+		                  min_measurements, push_res, 1);
 	write_reflist_2(tmp, split, sym);
 	reflist_free(split);
 
@@ -315,6 +316,7 @@ static void show_help(const char *s)
 "      --stop-after=<n>       Stop after merging <n> crystals.\n"
 "  -n, --iterations=<n>       Run <n> cycles of scaling and post-refinement.\n"
 "      --no-scale             Disable scale factor (G, B) refinement.\n"
+"      --no-Bscale            Disable B factor scaling.\n"
 "      --no-pr                Disable orientation/physics refinement.\n"
 "  -m, --model=<model>        Specify partiality model.\n"
 "      --min-measurements=<n> Minimum number of measurements to require.\n"
@@ -324,7 +326,12 @@ static void show_help(const char *s)
 "  -j <n>                     Run <n> analyses in parallel.\n"
 "      --no-free              Disable cross-validation (testing only).\n"
 "      --custom-split         List of files for custom dataset splitting.\n"
-"      --max-rel-B            Maximum allowable relative |B| factor.\n");
+"      --max-rel-B            Maximum allowable relative |B| factor.\n"
+"      --no-logs              Do not write extensive log files.\n"
+"  -w <pg>                    Apparent point group for resolving ambiguities.\n"
+"      --operator=<op>        Indexing ambiguity operator for resolving.\n"
+"      --force-bandwidth=<n>  Set all bandwidths to <n> (fraction).\n"
+"      --force-radius=<n>     Set all profile radii to <n> nm^-1.\n");
 }
 
 
@@ -578,7 +585,7 @@ static int set_initial_params(Crystal *cr, FILE *fh)
 
 	} else {
 
-		crystal_set_osf(cr, 0.0);
+		crystal_set_osf(cr, 1.0);
 		crystal_set_Bfac(cr, 0.0);
 
 	}
@@ -603,13 +610,22 @@ static void select_free_reflections(RefList *list, gsl_rng *rng)
 
 
 static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
-                            int fr)
+                            int fr, signed int inum)
 {
 	Reflection *refl;
 	RefListIterator *iter;
 	double G = crystal_get_osf(cr);
 	double B = crystal_get_Bfac(cr);
 	UnitCell *cell = crystal_get_cell(cr);
+	char ins[5];
+
+	if ( inum >= 0 ) {
+		snprintf(ins, 4, "%i", inum);
+	} else {
+		ins[0] = 'F';
+		ins[1] = '\0';
+	}
+
 
 	for ( refl = first_refl(list, &iter);
 	      refl != NULL;
@@ -622,6 +638,9 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
 
 		if ( !get_flag(refl) ) continue;  /* Not free-flagged */
 
+		/* Strong reflections only */
+		if ( get_intensity(refl) < 3.0*get_esd_intensity(refl) ) continue;
+
 		get_indices(refl, &h, &k, &l);
 		res = resolution(cell, h, k, l);
 		if ( 2.0*res > crystal_get_resolution_limit(cr) ) continue;
@@ -629,46 +648,55 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
 		match = find_refl(full, h, k, l);
 		if ( match == NULL ) continue;
 
+		/* Don't calculate pobs if reference reflection is weak */
+		if ( fabs(get_intensity(match)) / get_esd_intensity(match) < 3.0 ) continue;
+
 		/* Calculated partiality */
 		pcalc = get_partiality(refl);
 
 		/* Observed partiality */
-		corr = exp(-G) * exp(B*res*res) * get_lorentz(refl);
+		corr = G * exp(B*res*res) * get_lorentz(refl);
 		Ipart = get_intensity(refl) * corr;
 		pobs = Ipart / get_intensity(match);
 
-		fprintf(fh, "%5i %4i %4i %4i %8.4f %8.3f %8.3f\n",
-		        fr, h, k, l, 2*res/1e9, pcalc, pobs);
+		fprintf(fh, "%5i %4i %4i %4i %e %e %8.3f %8.3f %s\n",
+		        fr, h, k, l, 2*res, Ipart, pcalc, pobs, ins);
 
 	}
 }
 
 
 static void write_pgraph(RefList *full, Crystal **crystals, int n_crystals,
-                         int iter)
+                         signed int iter, const char *suff)
 {
 	FILE *fh;
 	char tmp[256];
 	int i;
 
-	snprintf(tmp, 256, "pgraph-iter%i.dat", iter);
+	snprintf(tmp, 256, "pr-logs/pgraph%s.dat", suff);
+
+	if ( iter == 0 ) {
+		fh = fopen(tmp, "w");
+	} else {
+		fh = fopen(tmp, "a");
+	}
 
-	fh = fopen(tmp, "w");
 	if ( fh == NULL ) {
 		ERROR("Failed to open '%s'\n", tmp);
 		return;
 	}
 
-	fprintf(fh, "   fr    h    k    l  1/d(nm)     pcalc   pobs\n");
+	if ( iter == 0 ) {
+		fprintf(fh, "  Crystal    h    k    l  1/d(m)   Ipart    pcalc   pobs   iteration\n");
+	}
 
 	for ( i=0; i<n_crystals; i++ ) {
 		if ( crystal_get_user_flag(crystals[i]) != 0 ) continue;
 		write_to_pgraph(fh, crystal_get_reflections(crystals[i]), full,
-		                crystals[i], i);
+		                crystals[i], i, iter);
 	}
 
 	fclose(fh);
-
 }
 
 
@@ -689,6 +717,7 @@ static void all_residuals(Crystal **crystals, int n_crystals, RefList *full,
 
 		if ( crystal_get_user_flag(crystals[i]) ) continue;
 
+		/* Scaling should have been done right before calling this */
 		r = residual(crystals[i], full, 0, NULL, NULL);
 		free_r = residual(crystals[i], full, 1, NULL, NULL);
 		log_r = log_residual(crystals[i], full, 0, NULL, NULL);
@@ -719,32 +748,82 @@ static void show_all_residuals(Crystal **crystals, int n_crystals,
 }
 
 
-static void dump_parameters(const char *filename, Crystal **crystals,
-                            int n_crystals)
+struct log_qargs
 {
-	FILE *fh;
-	fh = fopen(filename, "w");
-	if ( fh == NULL ) {
-		ERROR("Couldn't open partialator.params!\n");
-	} else {
-		fprintf(fh, "  cr        OSF       relB         div"
-		            " flag                      filename  event\n");
-		int i;
-		for ( i=0; i<n_crystals; i++ ) {
-			struct image *img;
-                        char *evt_str;
-			img = crystal_get_image(crystals[i]);
-                        evt_str = get_event_string(img->event);
-			fprintf(fh, "%4i %10.5f %10.2f %8.5e %-25s %s %s\n",
-			        i, crystal_get_osf(crystals[i]),
-			        crystal_get_Bfac(crystals[i])*1e20,
-			        crystal_get_image(crystals[i])->div,
-			        str_prflag(crystal_get_user_flag(crystals[i])),
-			        img->filename, evt_str);
-                        free(evt_str);
-		}
-		fclose(fh);
-	}
+	int iter;
+	int next;
+	Crystal **crystals;
+	int n_crystals;
+	RefList *full;
+	int scaleflags;
+};
+
+
+struct log_args
+{
+	Crystal *cr;
+	RefList *full;
+	int scaleflags;
+	int iter;
+	int cnum;
+};
+
+
+static void *get_log_task(void *vp)
+{
+	struct log_qargs *qargs = vp;
+	struct log_args *task;
+
+	if ( qargs->next >= qargs->n_crystals ) return NULL;
+
+	task = malloc(sizeof(struct log_args));
+	if ( task == NULL ) return NULL;
+
+	task->cr = qargs->crystals[qargs->next];
+	task->full = qargs->full;
+	task->iter = qargs->iter;
+	task->cnum = qargs->next;
+	task->scaleflags = qargs->scaleflags;
+
+	qargs->next += 20;
+	return task;
+}
+
+
+static void write_logs(void *vp, int cookie)
+{
+	struct log_args *args = vp;
+	write_specgraph(args->cr, args->full, args->iter, args->cnum);
+	write_gridscan(args->cr, args->full, args->iter, args->cnum,
+	               args->scaleflags);
+	write_test_logs(args->cr, args->full, args->iter, args->cnum);
+}
+
+
+static void done_log(void *qargs, void *vp)
+{
+	struct log_args *task = vp;
+	free(task);
+}
+
+
+static void write_logs_parallel(Crystal **crystals, int n_crystals,
+                                RefList *full, int iter, int n_threads,
+                                int scaleflags)
+{
+	struct log_qargs qargs;
+
+	qargs.iter = iter;
+	qargs.next = 0;
+	qargs.full = full;
+	qargs.crystals = crystals;
+	qargs.n_crystals = n_crystals;
+	qargs.scaleflags = scaleflags;
+
+	STATUS("Writing logs...\n");
+	run_threads(n_threads, write_logs, get_log_task, done_log, &qargs,
+	            n_crystals/20, 0, 0, 0);
+	STATUS("Done.\n");
 }
 
 
@@ -755,9 +834,10 @@ int main(int argc, char *argv[])
 	char *outfile = NULL;
 	char *sym_str = NULL;
 	SymOpList *sym;
+	SymOpList *amb;
+	SymOpList *w_sym;
 	int nthreads = 1;
 	int i;
-	struct image *images;
 	int n_iter = 10;
 	RefList *full;
 	int n_images = 0;
@@ -765,11 +845,12 @@ int main(int argc, char *argv[])
 	int n_crystals_seen = 0;
 	char cmdline[1024];
 	int no_scale = 0;
+	int no_Bscale = 0;
 	int no_pr = 0;
 	Stream *st;
 	Crystal **crystals;
 	char *pmodel_str = NULL;
-	PartialityModel pmodel = PMODEL_SCSPHERE;
+	PartialityModel pmodel = PMODEL_XSPHERE;
 	int min_measurements = 2;
 	char *rval;
 	int polarisation = 1;
@@ -785,6 +866,15 @@ int main(int argc, char *argv[])
 	char *csplit_fn = NULL;
 	struct custom_split *csplit = NULL;
 	double max_B = 1e-18;
+	char *rfile = NULL;
+	RefList *reference = NULL;
+	int no_logs = 0;
+	char *w_sym_str = NULL;
+	char *operator = NULL;
+	double force_bandwidth = -1.0;
+	double force_radius = -1.0;
+	char *audit_info;
+	int scaleflags = 0;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -805,10 +895,15 @@ int main(int argc, char *argv[])
 		{"push-res",           1, NULL,                5},
 		{"res-push",           1, NULL,                5}, /* compat */
 		{"custom-split",       1, NULL,                6},
-		{"max-rel-B"   ,       1, NULL,                7},
-		{"max-rel-b"   ,       1, NULL,                7}, /* compat */
+		{"max-rel-B",          1, NULL,                7},
+		{"max-rel-b",          1, NULL,                7}, /* compat */
+		{"reference",          1, NULL,                8}, /* ssshhh! */
+		{"operator",           1, NULL,                9},
+		{"force-bandwidth",    1, NULL,               10},
+		{"force-radius",       1, NULL,               11},
 
 		{"no-scale",           0, &no_scale,           1},
+		{"no-Bscale",          0, &no_Bscale,          1},
 		{"no-pr",              0, &no_pr,              1},
 		{"no-polarisation",    0, &polarisation,       0},
 		{"no-polarization",    0, &polarisation,       0}, /* compat */
@@ -816,6 +911,7 @@ int main(int argc, char *argv[])
 		{"polarization",       0, &polarisation,       1}, /* compat */
 		{"no-free",            0, &no_free,            1},
 		{"output-every-cycle", 0, &output_everycycle,  1},
+		{"no-logs",            0, &no_logs,            1},
 
 		{0, 0, NULL, 0}
 	};
@@ -827,7 +923,7 @@ int main(int argc, char *argv[])
 	}
 
 	/* Short options */
-	while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:j:m:",
+	while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:j:m:w:",
 	                        longopts, NULL)) != -1)
 	{
 
@@ -886,6 +982,10 @@ int main(int argc, char *argv[])
 			pmodel_str = strdup(optarg);
 			break;
 
+			case 'w' :
+			w_sym_str = strdup(optarg);
+			break;
+
 			case 2 :
 			errno = 0;
 			min_measurements = strtod(optarg, &rval);
@@ -932,6 +1032,33 @@ int main(int argc, char *argv[])
 			max_B = max_B * 1e-20;
 			break;
 
+			case 8 :
+			rfile = strdup(optarg);
+			break;
+
+			case 9 :
+			operator = strdup(optarg);
+			break;
+
+			case 10 :
+			errno = 0;
+			force_bandwidth = strtod(optarg, &rval);
+			if ( *rval != '\0' ) {
+				ERROR("Invalid value for --force-bandwidth.\n");
+				return 1;
+			}
+			break;
+
+			case 11 :
+			errno = 0;
+			force_radius = strtod(optarg, &rval);
+			if ( *rval != '\0' ) {
+				ERROR("Invalid value for --force-radius.\n");
+				return 1;
+			}
+			force_radius *= 1e9;
+			break;
+
 			case 0 :
 			break;
 
@@ -971,13 +1098,49 @@ int main(int argc, char *argv[])
 	free(sym_str);
 	if ( sym == NULL ) return 1;
 
+	if ( (w_sym_str != NULL) && (operator != NULL) ) {
+		ERROR("Specify the apparent symmetry (-w) or the operator, "
+		      "not both.\n");
+		return 1;
+	}
+
+	if ( w_sym_str == NULL ) {
+		w_sym = NULL;
+		amb = NULL;
+	} else {
+		pointgroup_warning(w_sym_str);
+		w_sym = get_pointgroup(w_sym_str);
+		free(w_sym_str);
+		if ( w_sym == NULL ) return 1;
+		amb = get_ambiguities(w_sym, sym);
+		if ( amb == NULL ) {
+			ERROR("Couldn't find ambiguity operator.\n");
+			ERROR("Check that your values for -y and -w are "
+			      "correct.\n");
+			return 1;
+		}
+
+	}
+
+	if ( operator ) {
+		amb = parse_symmetry_operations(operator);
+		if ( amb == NULL ) return 1;
+		set_symmetry_name(amb, "Ambiguity");
+	}
+
+	if ( amb != NULL ) {
+		STATUS("Indexing ambiguity resolution enabled.  "
+		       "The ambiguity operation(s) are:\n");
+		describe_symmetry(amb);
+		/* In contrast to ambigator, partialator can deal with multiple
+		 * ambiguities at once */
+	}
+
 	if ( pmodel_str != NULL ) {
 		if ( strcmp(pmodel_str, "unity") == 0 ) {
 			pmodel = PMODEL_UNITY;
-		} else if ( strcmp(pmodel_str, "scgaussian") == 0 ) {
-			pmodel = PMODEL_SCGAUSSIAN;
-		} else if ( strcmp(pmodel_str, "scsphere") == 0 ) {
-			pmodel = PMODEL_SCSPHERE;
+		} else if ( strcmp(pmodel_str, "xsphere") == 0 ) {
+			pmodel = PMODEL_XSPHERE;
 		} else if ( strcmp(pmodel_str, "random") == 0 ) {
 			pmodel = PMODEL_RANDOM;
 		} else {
@@ -988,7 +1151,25 @@ int main(int argc, char *argv[])
 
 	if ( (pmodel == PMODEL_UNITY) && !no_pr ) {
 		no_pr = 1;
-		STATUS("Setting --no-pr because we are not modelling partialities (--model=unity).\n");
+		STATUS("Setting --no-pr because we are not modelling "
+		       "partialities (--model=unity).\n");
+	}
+
+	if ( no_Bscale ) {
+		scaleflags |= SCALE_NO_B;
+	}
+
+	if ( !no_logs ) {
+		int r = mkdir("pr-logs", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+		if ( r ) {
+			if ( errno == EEXIST ) {
+				ERROR("WARNING: pr-logs folder already exists. "
+				      "Beware of mixing old and new log files!\n");
+			} else {
+				ERROR("Failed to create pr-logs folder.\n");
+				return 1;
+			}
+		}
 	}
 
 	/* Read the custom split list (if applicable) */
@@ -1001,6 +1182,20 @@ int main(int argc, char *argv[])
 		free(csplit_fn);
 	}
 
+	if ( rfile != NULL ) {
+		RefList *rread;
+		rread = read_reflections(rfile);
+		if ( rread == NULL ) {
+			ERROR("Failed to read reference reflections\n");
+			return 1;
+		}
+		reference = asymmetric_indices(rread, sym);
+		reflist_free(rread);
+		ERROR("WARNING: Using an external reference.\n");
+		ERROR("WARNING: If you publish a structure based on the result,"
+		      " expect to have to retract your paper!\n");
+	}
+
 	gsl_set_error_handler_off();
 	rng = gsl_rng_alloc(gsl_rng_mt19937);
 
@@ -1008,7 +1203,6 @@ int main(int argc, char *argv[])
 	n_images = 0;
 	n_crystals = 0;
 	n_crystals_seen = 0;
-	images = NULL;
 	crystals = NULL;
 	if ( sparams_fn != NULL ) {
 		char line[1024];
@@ -1032,35 +1226,27 @@ int main(int argc, char *argv[])
 
 		RefList *as;
 		int i;
-		struct image *images_new;
-		struct image *cur;
+		struct image cur;
 
-		images_new = realloc(images, (n_images+1)*sizeof(struct image));
-		if ( images_new == NULL ) {
-			ERROR("Failed to allocate memory for image list.\n");
-			return 1;
-		}
-		images = images_new;
-		cur = &images[n_images];
-
-		cur->div = NAN;
-		cur->bw = NAN;
-		cur->det = NULL;
-		if ( read_chunk_2(st, cur, STREAM_READ_REFLECTIONS
-		                           | STREAM_READ_UNITCELL) != 0 ) {
+		cur.div = NAN;
+		cur.bw = NAN;
+		cur.det = NULL;
+		if ( read_chunk_2(st, &cur, STREAM_READ_REFLECTIONS
+		                            | STREAM_READ_UNITCELL) != 0 ) {
 			break;
 		}
 
-		if ( isnan(cur->div) || isnan(cur->bw) ) {
+		if ( isnan(cur.div) || isnan(cur.bw) ) {
 			ERROR("Chunk doesn't contain beam parameters.\n");
 			return 1;
 		}
 
-		for ( i=0; i<cur->n_crystals; i++ ) {
+		for ( i=0; i<cur.n_crystals; i++ ) {
 
 			Crystal *cr;
 			Crystal **crystals_new;
 			RefList *cr_refl;
+			struct image *image;
 
 			n_crystals_seen++;
 			if ( n_crystals_seen <= start_after ) continue;
@@ -1073,11 +1259,19 @@ int main(int argc, char *argv[])
 				return 1;
 			}
 			crystals = crystals_new;
-			crystals[n_crystals] = cur->crystals[i];
+			crystals[n_crystals] = cur.crystals[i];
 			cr = crystals[n_crystals];
 
-			/* Image pointer will change due to later reallocs */
-			crystal_set_image(cr, NULL);
+			image = malloc(sizeof(struct image));
+			if ( image == NULL ) {
+				ERROR("Failed to allocatea memory for image.\n");
+				return 1;
+			}
+
+			crystal_set_image(cr, image);
+			*image = cur;
+			image->n_crystals = 1;
+			image->crystals = &crystals[n_crystals];
 
 			/* This is the raw list of reflections */
 			cr_refl = crystal_get_reflections(cr);
@@ -1087,7 +1281,7 @@ int main(int argc, char *argv[])
 			if ( polarisation ) {
 				polarisation_correction(cr_refl,
 						        crystal_get_cell(cr),
-						        cur);
+						        image);
 			}
 
 			if ( !no_free ) select_free_reflections(cr_refl, rng);
@@ -1120,68 +1314,77 @@ int main(int argc, char *argv[])
 	display_progress(n_images, n_crystals);
 	fprintf(stderr, "\n");
 	if ( sparams_fh != NULL ) fclose(sparams_fh);
-
+	audit_info = stream_audit_info(st);
 	close_stream(st);
 
-	/* Fill in image pointers */
-	for ( i=0; i<n_images; i++ ) {
-		int j;
-		for ( j=0; j<images[i].n_crystals; j++ ) {
-
-			Crystal *cryst;
-
-			cryst = images[i].crystals[j];
-			crystal_set_image(cryst, &images[i]);
-
-			/* Now it's safe to do the following */
-			update_predictions(cryst);
-			calculate_partialities(cryst, pmodel);
-
+	STATUS("Initial partiality calculation...\n");
+	for ( i=0; i<n_crystals; i++ ) {
+		Crystal *cr = crystals[i];
+		if ( force_bandwidth > 0.0 ) {
+			crystal_get_image(cr)->bw = force_bandwidth;
 		}
+		if ( force_radius > 0.0 ) {
+			crystal_set_profile_radius(cr, force_radius);
+		}
+		update_predictions(cr);
+		calculate_partialities(cr, pmodel);
 	}
 
 	if (csplit != NULL) check_csplit(crystals, n_crystals, csplit);
 
 	/* Make a first pass at cutting out crap */
-	STATUS("Checking patterns.\n");
+	//STATUS("Early rejection...\n");
 	//early_rejection(crystals, n_crystals);
 
-	/* Initial rejection, figures of merit etc */
-	full = merge_intensities(crystals, n_crystals, nthreads, pmodel,
-	                         min_measurements, push_res, 1);
+	/* Create reference data set if we don't already have one */
+	if ( reference == NULL ) {
+		if ( !no_scale ) {
+			STATUS("Initial scaling...\n");
+			scale_all(crystals, n_crystals, nthreads, scaleflags);
+		}
+		full = merge_intensities(crystals, n_crystals, nthreads,
+		                         min_measurements, push_res, 1);
+	} else {
+		full = reference;
+	}
+
+	/* Check rejection and write figures of merit */
 	check_rejection(crystals, n_crystals, full, max_B);
 	show_all_residuals(crystals, n_crystals, full);
-	write_pgraph(full, crystals, n_crystals, 0);
+	if ( !no_pr && !no_logs ) {
+		write_pgraph(full, crystals, n_crystals, 0, "");
+		write_logs_parallel(crystals, n_crystals, full, 0, nthreads,
+		                    scaleflags);
+	}
 
 	/* Iterate */
 	for ( i=0; i<n_iter; i++ ) {
 
 		STATUS("Scaling and refinement cycle %i of %i\n", i+1, n_iter);
 
-		if ( !no_scale ) {
-			scale_all(crystals, n_crystals, nthreads, pmodel);
-			reflist_free(full);
-			full = merge_intensities(crystals, n_crystals, nthreads,
-			                         pmodel, min_measurements,
-			                         push_res, 1);
+		if ( !no_pr ) {
+			refine_all(crystals, n_crystals, full, nthreads, pmodel,
+			           0, i+1, no_logs, sym, amb, scaleflags);
 		}
 
-		if ( !no_pr ) {
-			refine_all(crystals, n_crystals, full, nthreads,
-			           pmodel);
+		/* Create new reference if needed */
+		if ( reference == NULL ) {
 			reflist_free(full);
+			if ( !no_scale ) {
+				scale_all(crystals, n_crystals, nthreads,
+				          scaleflags);
+			}
 			full = merge_intensities(crystals, n_crystals, nthreads,
-			                         pmodel, min_measurements,
+			                         min_measurements,
 			                         push_res, 1);
-		}
+		} /* else full still equals reference */
 
 		check_rejection(crystals, n_crystals, full, max_B);
-		reflist_free(full);
-		full = merge_intensities(crystals, n_crystals, nthreads,
-		                         pmodel, min_measurements,
-		                         push_res, 1);
 		show_all_residuals(crystals, n_crystals, full);
-		write_pgraph(full, crystals, n_crystals, i+1);
+
+		if ( !no_logs ) {
+			write_pgraph(full, crystals, n_crystals, i+1, "");
+		}
 
 		if ( output_everycycle ) {
 
@@ -1208,21 +1411,37 @@ int main(int argc, char *argv[])
 				}
 			}
 
-			/* Dump parameters */
-			snprintf(tmp, 1024, "iter%.2d_partialator.params", i+1);
-			dump_parameters(tmp, crystals, n_crystals);
+		}
+	}
 
+	/* Final merge */
+	STATUS("Final merge...\n");
+	if ( reference == NULL ) {
+		reflist_free(full);
+		if ( !no_scale ) {
+			scale_all(crystals, n_crystals, nthreads, scaleflags);
 		}
+		full = merge_intensities(crystals, n_crystals, nthreads,
+		                         min_measurements,
+		                         push_res, 1);
+	} else {
+		full = merge_intensities(crystals, n_crystals, nthreads,
+		                         min_measurements, push_res, 1);
 	}
 
-	full = merge_intensities(crystals, n_crystals, nthreads, pmodel,
-	                         min_measurements, push_res, 1);
+	/* Write final figures of merit (no rejection any more) */
 	show_all_residuals(crystals, n_crystals, full);
-	write_pgraph(full, crystals, n_crystals, n_iter+1);
+	if ( !no_pr && !no_logs ) {
+		write_pgraph(full, crystals, n_crystals, -1, "");
+		write_logs_parallel(crystals, n_crystals, full, -1, nthreads,
+		                    scaleflags);
+	}
 
 	/* Output results */
 	STATUS("Writing overall results to %s\n", outfile);
 	reflist_add_command_and_version(full, argc, argv);
+	reflist_add_notes(full, "Audit information from stream:");
+	reflist_add_notes(full, audit_info);
 	write_reflist_2(outfile, full, sym);
 
 	/* Output split results */
@@ -1238,9 +1457,6 @@ int main(int argc, char *argv[])
 		}
 	}
 
-	/* Dump parameters */
-	dump_parameters("partialator.params", crystals, n_crystals);
-
 	/* Clean up */
 	gsl_rng_free(rng);
 	for ( i=0; i<n_crystals; i++ ) {
@@ -1248,13 +1464,9 @@ int main(int argc, char *argv[])
 		crystal_free(crystals[i]);
 	}
 	reflist_free(full);
-	free(sym);
+	free_symoplist(sym);
 	free(outfile);
 	free(crystals);
-	for ( i=0; i<n_images; i++ ) {
-		free(images[i].filename);
-	}
-	free(images);
 	free(infile);
 
 	return 0;