partialator.c 12.7 KB
Newer Older
1
/*
2
 * partialator.c
3
 *
4
 * Scaling and post refinement for coherent nanocrystallography
5
 *
Thomas White's avatar
Thomas White committed
6
 * (c) 2006-2011 Thomas White <taw@physics.org>
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
 *
 * Part of CrystFEL - crystallography with a FEL
 *
 */


#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdarg.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <unistd.h>
#include <getopt.h>
Thomas White's avatar
Thomas White committed
23
#include <assert.h>
Thomas White's avatar
Tidy up    
Thomas White committed
24
#include <pthread.h>
25
#include <gsl/gsl_errno.h>
26

Thomas White's avatar
Thomas White committed
27
#include "utils.h"
28
#include "hdf5-file.h"
Thomas White's avatar
Thomas White committed
29
#include "symmetry.h"
Thomas White's avatar
Thomas White committed
30
#include "stream.h"
31
#include "geometry.h"
32
#include "peaks.h"
Thomas White's avatar
Thomas White committed
33
#include "thread-pool.h"
Thomas White's avatar
Thomas White committed
34
#include "beam-parameters.h"
35
#include "post-refinement.h"
36
#include "hrs-scaling.h"
Thomas White's avatar
Thomas White committed
37
#include "reflist.h"
Thomas White's avatar
Thomas White committed
38
#include "reflist-utils.h"
Thomas White's avatar
Thomas White committed
39
#include "scaling-report.h"
40
41
42
43
44
45


static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
46
"Scaling and post refinement for coherent nanocrystallography.\n"
47
48
49
"\n"
"  -h, --help                 Display this help message.\n"
"\n"
Thomas White's avatar
Thomas White committed
50
51
"  -i, --input=<filename>     Specify the name of the input 'stream'.\n"
"                              (must be a file, not e.g. stdin)\n"
Thomas White's avatar
Thomas White committed
52
"  -o, --output=<filename>    Output filename.  Default: facetron.hkl.\n"
Thomas White's avatar
Thomas White committed
53
"  -g. --geometry=<file>      Get detector geometry from file.\n"
54
55
56
57
"  -b, --beam=<file>          Get beam parameters from file, which provides\n"
"                              initial values for parameters, and nominal\n"
"                              wavelengths if no per-shot value is found in \n"
"                              an HDF5 file.\n"
58
"  -y, --symmetry=<sym>       Merge according to symmetry <sym>.\n"
59
"  -n, --iterations=<n>       Run <n> cycles of scaling and post-refinement.\n"
60
61
62
"      --reference=<file>     Refine images against reflections in <file>,\n"
"                              instead of taking the mean of the intensity\n"
"                              estimates.\n"
63
"\n"
Thomas White's avatar
Thomas White committed
64
65
66
67
"  -j <n>                     Run <n> analyses in parallel.\n");
}


Thomas White's avatar
Thomas White committed
68
struct refine_args
Thomas White's avatar
Thomas White committed
69
{
70
	RefList *full;
Thomas White's avatar
Thomas White committed
71
	struct image *image;
72
	FILE *graph;
73
	FILE *pgraph;
Thomas White's avatar
Thomas White committed
74
75
76
};


Thomas White's avatar
Thomas White committed
77
struct queue_args
Thomas White's avatar
Thomas White committed
78
{
Thomas White's avatar
Thomas White committed
79
80
81
82
83
84
85
86
87
88
89
	int n;
	int n_done;
	int n_total_patterns;
	struct image *images;
	struct refine_args task_defaults;
};


static void refine_image(void *task, int id)
{
	struct refine_args *pargs = task;
90
	struct image *image = pargs->image;
Thomas White's avatar
Thomas White committed
91
	image->id = id;
92

93
	pr_refine(image, pargs->full);
94
95
96
}


Thomas White's avatar
Thomas White committed
97
static void *get_image(void *vqargs)
98
{
Thomas White's avatar
Thomas White committed
99
100
	struct refine_args *task;
	struct queue_args *qargs = vqargs;
Thomas White's avatar
Thomas White committed
101

Thomas White's avatar
Thomas White committed
102
103
	task = malloc(sizeof(struct refine_args));
	memcpy(task, &qargs->task_defaults, sizeof(struct refine_args));
104

Thomas White's avatar
Thomas White committed
105
	task->image = &qargs->images[qargs->n];
106

Thomas White's avatar
Thomas White committed
107
108
109
110
	qargs->n++;

	return task;
}
Thomas White's avatar
Thomas White committed
111

112

Thomas White's avatar
Thomas White committed
113
114
115
116
117
118
static void done_image(void *vqargs, void *task)
{
	struct queue_args *qargs = vqargs;

	qargs->n_done++;

Thomas White's avatar
Thomas White committed
119
	progress_bar(qargs->n_done, qargs->n_total_patterns, "Refining");
Thomas White's avatar
Thomas White committed
120
121
122
123
124
	free(task);
}


static void refine_all(struct image *images, int n_total_patterns,
125
126
                       struct detector *det,
                       RefList *full, int nthreads,
Thomas White's avatar
Thomas White committed
127
128
129
130
131
132
133
134
135
136
137
138
139
                       FILE *graph, FILE *pgraph)
{
	struct refine_args task_defaults;
	struct queue_args qargs;

	task_defaults.full = full;
	task_defaults.image = NULL;
	task_defaults.graph = graph;
	task_defaults.pgraph = pgraph;

	qargs.task_defaults = task_defaults;
	qargs.n = 0;
	qargs.n_done = 0;
Thomas White's avatar
Thomas White committed
140
141
	qargs.n_total_patterns = n_total_patterns;
	qargs.images = images;
Thomas White's avatar
Thomas White committed
142

Thomas White's avatar
Thomas White committed
143
	/* Don't have threads which are doing nothing */
Thomas White's avatar
Whoops    
Thomas White committed
144
	if ( n_total_patterns < nthreads ) nthreads = n_total_patterns;
Thomas White's avatar
Thomas White committed
145

Thomas White's avatar
Thomas White committed
146
	run_threads(nthreads, refine_image, get_image, done_image,
Thomas White's avatar
Thomas White committed
147
	            &qargs, n_total_patterns, 0, 0, 0);
148
149
150
}


151
/* Decide which reflections can be scaled */
152
static int select_scalable_reflections(RefList *list, RefList *reference)
153
154
155
156
157
158
159
160
161
{
	Reflection *refl;
	RefListIterator *iter;
	int nobs = 0;

	for ( refl = first_refl(list, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {

162
		int sc = 1;
163
164
		double v;

165
166
167
		/* This means the reflection was not found on the last check */
		if ( get_redundancy(refl) == 0 ) sc = 0;

168
		if ( get_partiality(refl) < 0.1 ) sc = 0;
169
		v = fabs(get_intensity(refl));
170
		if ( v < 0.1 ) sc = 0;
171

172
173
174
		/* If we are scaling against a reference set, we additionally
		 * require that this reflection is in the reference list. */
		if ( reference != NULL ) {
175
176
			signed int h, k, l;
			get_indices(refl, &h, &k, &l);
177
			if ( find_refl(reference, h, k, l) == NULL ) sc = 0;
178
179
		}

180
181
182
		set_scalable(refl, sc);

		if ( sc ) nobs++;
183
184
185
186
187
188
	}

	return nobs;
}


189
static void select_reflections_for_refinement(struct image *images, int n,
190
                                              RefList *full, int have_reference)
191
192
193
194
195
196
197
{
	int i;

	for ( i=0; i<n; i++ ) {

		Reflection *refl;
		RefListIterator *iter;
198
199
200
201
202
		int n_acc = 0;
		int n_nomatch = 0;
		int n_noscale = 0;
		int n_fewmatch = 0;
		int n_ref = 0;
203
204
205
206
207
208
209
210

		for ( refl = first_refl(images[i].reflections, &iter);
		      refl != NULL;
		      refl = next_refl(refl, iter) )
		{
			signed int h, k, l;
			int sc;

211
212
213
214
215
216
			n_ref++;

			/* We require that the reflection itself is scalable
			 * (i.e. sensible partiality and intensity) and that
			 * the "full" estimate of this reflection is made from
			 * at least two parts. */
217
218
			get_indices(refl, &h, &k, &l);
			sc = get_scalable(refl);
219
220
221
			if ( !sc ) {

				n_noscale++;
222
				set_refinable(refl, 0);
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239

			} else {

				Reflection *f = find_refl(full, h, k, l);

				if ( f != NULL ) {

					int r = get_redundancy(f);
					if ( (r >= 2) || have_reference ) {
						set_refinable(refl, 1);
						n_acc++;
					} else {
						n_fewmatch++;
					}

				} else {
					n_nomatch++;
240
					set_refinable(refl, 0);
241
				}
242
243
244
245

			}
		}

246
247
248
249
250
251
252
253
		STATUS("Image %4i: %i guide reflections accepted "
		       "(%i not scalable, %i few matches, %i total)\n",
		       i, n_acc, n_noscale, n_fewmatch, n_ref);

		/* This would be a silly situation, since there must be a match
		 * if THIS pattern has a scalable part of the reflection! */
		assert(n_nomatch == 0);

254
255
256
257
	}
}


Thomas White's avatar
Thomas White committed
258
259
260
261
int main(int argc, char *argv[])
{
	int c;
	char *infile = NULL;
Thomas White's avatar
Thomas White committed
262
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
263
	char *geomfile = NULL;
Thomas White's avatar
Thomas White committed
264
265
	char *sym = NULL;
	FILE *fh;
Thomas White's avatar
Thomas White committed
266
267
	int nthreads = 1;
	struct detector *det;
Thomas White's avatar
Thomas White committed
268
269
	int i;
	int n_total_patterns;
270
271
	struct image *images;
	int n_iter = 10;
Thomas White's avatar
Thomas White committed
272
	struct beam_params *beam = NULL;
273
	RefList *full;
274
275
276
	int n_found = 0;
	int n_expected = 0;
	int n_notfound = 0;
Thomas White's avatar
Thomas White committed
277
	int n_usable_patterns = 0;
278
	int nobs;
279
	char *reference_file = NULL;
280
	RefList *reference = NULL;
281
	int n_dud;
282
	int have_reference = 0;
Thomas White's avatar
Thomas White committed
283
284
285
286
287

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
288
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
289
		{"geometry",           1, NULL,               'g'},
Thomas White's avatar
Thomas White committed
290
		{"beam",               1, NULL,               'b'},
291
		{"symmetry",           1, NULL,               'y'},
292
		{"iterations",         1, NULL,               'n'},
293
		{"reference",          1, NULL,                1},
Thomas White's avatar
Thomas White committed
294
295
296
297
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
298
	while ((c = getopt_long(argc, argv, "hi:g:x:j:y:o:b:",
299
300
	                        longopts, NULL)) != -1)
	{
Thomas White's avatar
Thomas White committed
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318

		switch (c) {
		case 'h' :
			show_help(argv[0]);
			return 0;

		case 'i' :
			infile = strdup(optarg);
			break;

		case 'g' :
			geomfile = strdup(optarg);
			break;

		case 'j' :
			nthreads = atoi(optarg);
			break;

319
320
321
322
		case 'y' :
			sym = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
323
324
325
326
		case 'o' :
			outfile = strdup(optarg);
			break;

327
328
329
330
		case 'n' :
			n_iter = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
331
332
333
334
335
336
337
338
339
		case 'b' :
			beam = get_beam_parameters(optarg);
			if ( beam == NULL ) {
				ERROR("Failed to load beam parameters"
				      " from '%s'\n", optarg);
				return 1;
			}
			break;

340
		case 1 :
341
			reference_file = strdup(optarg);
342
343
			break;

Thomas White's avatar
Thomas White committed
344
345
346
347
348
349
350
351
352
		case 0 :
			break;

		default :
			return 1;
		}

	}

Thomas White's avatar
Thomas White committed
353
	/* Sanitise input filename and open */
Thomas White's avatar
Thomas White committed
354
355
356
357
358
359
360
361
362
363
364
365
366
	if ( infile == NULL ) {
		infile = strdup("-");
	}
	if ( strcmp(infile, "-") == 0 ) {
		fh = stdin;
	} else {
		fh = fopen(infile, "r");
	}
	if ( fh == NULL ) {
		ERROR("Failed to open input file '%s'\n", infile);
		return 1;
	}

Thomas White's avatar
Thomas White committed
367
368
	/* Sanitise output filename */
	if ( outfile == NULL ) {
Thomas White's avatar
Whoops    
Thomas White committed
369
		outfile = strdup("partialator.hkl");
Thomas White's avatar
Thomas White committed
370
371
	}

Thomas White's avatar
Thomas White committed
372
373
	if ( sym == NULL ) sym = strdup("1");

Thomas White's avatar
Thomas White committed
374
	/* Get detector geometry */
Thomas White's avatar
Thomas White committed
375
376
377
378
379
380
381
	det = get_detector_geometry(geomfile);
	if ( det == NULL ) {
		ERROR("Failed to read detector geometry from '%s'\n", geomfile);
		return 1;
	}
	free(geomfile);

Thomas White's avatar
Thomas White committed
382
383
384
385
386
	if ( beam == NULL ) {
		ERROR("You must provide a beam parameters file.\n");
		return 1;
	}

387
	if ( reference_file != NULL ) {
Thomas White's avatar
Thomas White committed
388

389
		RefList *list;
390

391
		list = read_reflections(reference_file);
392
		free(reference_file);
393
		if ( list == NULL ) return 1;
394
		reference = asymmetric_indices(list, sym);
395
		reflist_free(list);
396
		have_reference = 1;
397

398
399
	}

Thomas White's avatar
Thomas White committed
400
	n_total_patterns = count_patterns(fh);
Thomas White's avatar
Thomas White committed
401
402
403
404
	if ( n_total_patterns == 0 ) {
		ERROR("No patterns to process.\n");
		return 1;
	}
Thomas White's avatar
Thomas White committed
405
406
	STATUS("There are %i patterns to process\n", n_total_patterns);

407
408
	gsl_set_error_handler_off();

409
410
411
412
413
414
415
416
	images = malloc(n_total_patterns * sizeof(struct image));
	if ( images == NULL ) {
		ERROR("Couldn't allocate memory for images.\n");
		return 1;
	}

	/* Fill in what we know about the images so far */
	rewind(fh);
417
	nobs = 0;
418
419
	for ( i=0; i<n_total_patterns; i++ ) {

420
		RefList *as;
Thomas White's avatar
Thomas White committed
421
		struct image *cur;
422
		int nn_expected, nn_found, nn_notfound;
423

Thomas White's avatar
Thomas White committed
424
		cur = &images[n_usable_patterns];
425

Thomas White's avatar
Thomas White committed
426
427
428
		cur->det = det;

		if ( read_chunk(fh, cur) != 0 ) {
Thomas White's avatar
Thomas White committed
429
430
			/* Should not happen, because we counted the patterns
			 * earlier. */
Thomas White's avatar
Thomas White committed
431
			ERROR("Failed to read chunk from the input stream.\n");
432
433
434
			return 1;
		}

435
		/* Won't be needing this, if it exists */
Thomas White's avatar
Thomas White committed
436
437
		image_feature_list_free(cur->features);
		cur->features = NULL;
438

Thomas White's avatar
Thomas White committed
439
		/* "n_usable_patterns" will not be incremented in this case */
Thomas White's avatar
Thomas White committed
440
		if ( cur->indexed_cell == NULL ) continue;
441

Thomas White's avatar
Thomas White committed
442
		/* Fill in initial estimates of stuff */
Thomas White's avatar
Thomas White committed
443
444
445
446
447
		cur->div = beam->divergence;
		cur->bw = beam->bandwidth;
		cur->width = det->max_fs;
		cur->height = det->max_ss;
		cur->osf = 1.0;
Thomas White's avatar
Thomas White committed
448
		cur->profile_radius = 0.003e9;
Thomas White's avatar
Thomas White committed
449
		cur->pr_dud = 0;
450

Thomas White's avatar
Thomas White committed
451
		/* Muppet proofing */
Thomas White's avatar
Thomas White committed
452
453
454
		cur->data = NULL;
		cur->flags = NULL;
		cur->beam = NULL;
455

Thomas White's avatar
Thomas White committed
456
		/* This is the raw list of reflections */
Thomas White's avatar
Thomas White committed
457
458
459
		as = asymmetric_indices(cur->reflections, sym);
		reflist_free(cur->reflections);
		cur->reflections = as;
460

461
462
463
464
		update_partialities(cur, &nn_expected, &nn_found, &nn_notfound);
		n_expected += nn_expected;
		n_found += nn_found;
		n_notfound += nn_notfound;
465

466
467
		nobs += select_scalable_reflections(cur->reflections,
		                                    reference);
468

469
		progress_bar(i, n_total_patterns-1, "Loading pattern data");
Thomas White's avatar
Thomas White committed
470
		n_usable_patterns++;
471
472
473

	}
	fclose(fh);
Thomas White's avatar
Thomas White committed
474
475
	STATUS("Found %5.2f%% of the expected peaks (missed %i of %i).\n",
	       100.0 * (double)n_found / n_expected, n_notfound, n_expected);
Thomas White's avatar
Thomas White committed
476

477
	/* Make initial estimates */
478
	STATUS("Performing initial scaling.\n");
479
	full = scale_intensities(images, n_usable_patterns, reference);
Thomas White's avatar
Thomas White committed
480

Thomas White's avatar
Thomas White committed
481
	/* Iterate */
482
	for ( i=0; i<n_iter; i++ ) {
Thomas White's avatar
Thomas White committed
483

484
485
		FILE *fhg;
		FILE *fhp;
486
		char filename[1024];
487
		int j;
488
		RefList *comp;
489

490
		STATUS("Post refinement cycle %i of %i\n", i+1, n_iter);
Thomas White's avatar
Thomas White committed
491

492
493
494
495
496
497
498
499
500
501
		snprintf(filename, 1023, "p-iteration-%i.dat", i+1);
		fhg = fopen(filename, "w");
		if ( fhg == NULL ) {
			ERROR("Failed to open '%s'\n", filename);
			/* Nothing will be written later */
		}

		snprintf(filename, 1023, "g-iteration-%i.dat", i+1);
		fhp = fopen(filename, "w");
		if ( fhp == NULL ) {
502
503
504
505
			ERROR("Failed to open '%s'\n", filename);
			/* Nothing will be written later */
		}

506
507
508
509
510
		if ( reference == NULL ) {
			comp = full;
		} else {
			comp = reference;
		}
511

Thomas White's avatar
Thomas White committed
512
		/* Refine the geometry of all patterns to get the best fit */
513
514
		select_reflections_for_refinement(images, n_usable_patterns,
		                                  comp, have_reference);
515
516
		refine_all(images, n_usable_patterns, det,
		           comp, nthreads, fhg, fhp);
Thomas White's avatar
Thomas White committed
517

518
		nobs = 0;
519
520
521
522
		for ( j=0; j<n_usable_patterns; j++ ) {

			struct image *cur = &images[j];

523
			nobs += select_scalable_reflections(cur->reflections,
524
			                                    reference);
525
526
527

		}

Thomas White's avatar
Thomas White committed
528
		/* Re-estimate all the full intensities */
529
		reflist_free(full);
Thomas White's avatar
Thomas White committed
530
		full = scale_intensities(images, n_usable_patterns,
531
532
		                         reference);

533
534
		fclose(fhg);
		fclose(fhp);
Thomas White's avatar
Thomas White committed
535
536
537
538

	}

	STATUS("Final scale factors:\n");
539
	n_dud = 0;
Thomas White's avatar
Thomas White committed
540
	for ( i=0; i<n_usable_patterns; i++ ) {
541
		if ( images[i].pr_dud ) n_dud++;
Thomas White's avatar
Thomas White committed
542
		STATUS("%4i : %5.2f\n", i, images[i].osf);
Thomas White's avatar
Thomas White committed
543
	}
544
	STATUS("%i images could not be refined on the last cycle.\n", n_dud);
Thomas White's avatar
Thomas White committed
545
546

	/* Output results */
547
	write_reflist(outfile, full, images[0].indexed_cell);
Thomas White's avatar
Thomas White committed
548

Thomas White's avatar
Thomas White committed
549
550
	scaling_report("scaling-report.pdf", images, n_usable_patterns, infile);

Thomas White's avatar
Thomas White committed
551
	/* Clean up */
Thomas White's avatar
Thomas White committed
552
	for ( i=0; i<n_usable_patterns; i++ ) {
Thomas White's avatar
Thomas White committed
553
		reflist_free(images[i].reflections);
Thomas White's avatar
Thomas White committed
554
	}
555
	reflist_free(full);
Thomas White's avatar
Thomas White committed
556
	free(sym);
Thomas White's avatar
Thomas White committed
557
	free(outfile);
Thomas White's avatar
Thomas White committed
558
	free_detector_geometry(det);
559
	free(beam);
560
	if ( reference != NULL ) {
561
		reflist_free(reference);
562
	}
Thomas White's avatar
Thomas White committed
563
	for ( i=0; i<n_usable_patterns; i++ ) {
564
565
566
567
		cell_free(images[i].indexed_cell);
		free(images[i].filename);
	}
	free(images);
Thomas White's avatar
Thomas White committed
568
	free(infile);
569
570
571

	return 0;
}