partialator.c 12.9 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"
39
40
41
42
43
44


static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
45
"Scaling and post refinement for coherent nanocrystallography.\n"
46
47
48
"\n"
"  -h, --help                 Display this help message.\n"
"\n"
Thomas White's avatar
Thomas White committed
49
50
"  -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
51
"  -o, --output=<filename>    Output filename.  Default: facetron.hkl.\n"
Thomas White's avatar
Thomas White committed
52
"  -g. --geometry=<file>      Get detector geometry from file.\n"
53
54
55
56
"  -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"
57
"  -y, --symmetry=<sym>       Merge according to symmetry <sym>.\n"
58
"  -n, --iterations=<n>       Run <n> cycles of scaling and post-refinement.\n"
59
60
61
"      --reference=<file>     Refine images against reflections in <file>,\n"
"                              instead of taking the mean of the intensity\n"
"                              estimates.\n"
62
"\n"
Thomas White's avatar
Thomas White committed
63
64
65
66
"  -j <n>                     Run <n> analyses in parallel.\n");
}


Thomas White's avatar
Thomas White committed
67
struct refine_args
Thomas White's avatar
Thomas White committed
68
{
Thomas White's avatar
Thomas White committed
69
70
	const char *sym;
	ReflItemList *obs;
71
	RefList *full;
Thomas White's avatar
Thomas White committed
72
	struct image *image;
73
	FILE *graph;
74
	FILE *pgraph;
Thomas White's avatar
Thomas White committed
75
76
77
};


Thomas White's avatar
Thomas White committed
78
struct queue_args
Thomas White's avatar
Thomas White committed
79
{
Thomas White's avatar
Thomas White committed
80
81
82
83
84
85
86
87
88
89
90
	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;
91
	struct image *image = pargs->image;
Thomas White's avatar
Thomas White committed
92
	image->id = id;
93

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


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

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

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

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

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

113

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

	qargs->n_done++;

Thomas White's avatar
Thomas White committed
120
	progress_bar(qargs->n_done, qargs->n_total_patterns, "Refining");
Thomas White's avatar
Thomas White committed
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
	free(task);
}


static void refine_all(struct image *images, int n_total_patterns,
                       struct detector *det, const char *sym,
                       ReflItemList *obs, RefList *full, int nthreads,
                       FILE *graph, FILE *pgraph)
{
	struct refine_args task_defaults;
	struct queue_args qargs;

	task_defaults.sym = sym;
	task_defaults.obs = obs;
	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
143
144
	qargs.n_total_patterns = n_total_patterns;
	qargs.images = images;
Thomas White's avatar
Thomas White committed
145

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

Thomas White's avatar
Thomas White committed
149
	run_threads(nthreads, refine_image, get_image, done_image,
Thomas White's avatar
Thomas White committed
150
	            &qargs, n_total_patterns, 0, 0, 0);
151
152
153
}


154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
/* Decide which reflections can be scaled */
static int select_scalable_reflections(RefList *list, ReflItemList *sc_l)
{
	Reflection *refl;
	RefListIterator *iter;
	int nobs = 0;

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

		int scalable = 1;
		double v;

		if ( get_partiality(refl) < 0.1 ) scalable = 0;
		v = fabs(get_intensity(refl));
		if ( v < 0.1 ) scalable = 0;
		set_scalable(refl, scalable);

		if ( scalable ) {

			signed int h, k, l;

			nobs++;

			/* Add (asymmetric) indices to list */
			get_indices(refl, &h, &k, &l);

			if ( !find_item(sc_l, h, k, l) ) {
				add_item(sc_l, h, k, l);
			}

		}

	}

	return nobs;
}


194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
static void select_reflections_for_refinement(struct image *images, int n,
                                              ReflItemList *scalable)
{
	int i;

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

		Reflection *refl;
		RefListIterator *iter;

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

			get_indices(refl, &h, &k, &l);
			sc = get_scalable(refl);

			if ( sc && find_item(scalable, h, k, l) ) {
				set_refinable(refl, 1);
			}
		}

	}
}


Thomas White's avatar
Thomas White committed
223
224
225
226
int main(int argc, char *argv[])
{
	int c;
	char *infile = NULL;
Thomas White's avatar
Thomas White committed
227
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
228
	char *geomfile = NULL;
Thomas White's avatar
Thomas White committed
229
230
	char *sym = NULL;
	FILE *fh;
Thomas White's avatar
Thomas White committed
231
232
	int nthreads = 1;
	struct detector *det;
233
	ReflItemList *scalable;
Thomas White's avatar
Thomas White committed
234
235
	int i;
	int n_total_patterns;
236
237
	struct image *images;
	int n_iter = 10;
Thomas White's avatar
Thomas White committed
238
	struct beam_params *beam = NULL;
239
	RefList *full;
240
241
242
	int n_found = 0;
	int n_expected = 0;
	int n_notfound = 0;
243
	char *cref;
Thomas White's avatar
Thomas White committed
244
	int n_usable_patterns = 0;
245
	int nobs;
246
	char *reference_file = NULL;
247
	double *reference = NULL;
248
	RefList *reference_list = NULL;
249
	int n_dud;
Thomas White's avatar
Thomas White committed
250
251
252
253
254

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
255
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
256
		{"geometry",           1, NULL,               'g'},
Thomas White's avatar
Thomas White committed
257
		{"beam",               1, NULL,               'b'},
258
		{"symmetry",           1, NULL,               'y'},
259
		{"iterations",         1, NULL,               'n'},
260
		{"reference",          1, NULL,                1},
Thomas White's avatar
Thomas White committed
261
262
263
264
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
265
	while ((c = getopt_long(argc, argv, "hi:g:x:j:y:o:b:",
266
267
	                        longopts, NULL)) != -1)
	{
Thomas White's avatar
Thomas White committed
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285

		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;

286
287
288
289
		case 'y' :
			sym = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
290
291
292
293
		case 'o' :
			outfile = strdup(optarg);
			break;

294
295
296
297
		case 'n' :
			n_iter = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
298
299
300
301
302
303
304
305
306
		case 'b' :
			beam = get_beam_parameters(optarg);
			if ( beam == NULL ) {
				ERROR("Failed to load beam parameters"
				      " from '%s'\n", optarg);
				return 1;
			}
			break;

307
		case 1 :
308
			reference_file = strdup(optarg);
309
310
			break;

Thomas White's avatar
Thomas White committed
311
312
313
314
315
316
317
318
319
		case 0 :
			break;

		default :
			return 1;
		}

	}

Thomas White's avatar
Thomas White committed
320
	/* Sanitise input filename and open */
Thomas White's avatar
Thomas White committed
321
322
323
324
325
326
327
328
329
330
331
332
333
334
	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;
	}
	free(infile);

Thomas White's avatar
Thomas White committed
335
336
	/* Sanitise output filename */
	if ( outfile == NULL ) {
Thomas White's avatar
Whoops    
Thomas White committed
337
		outfile = strdup("partialator.hkl");
Thomas White's avatar
Thomas White committed
338
339
	}

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

Thomas White's avatar
Thomas White committed
342
	/* Get detector geometry */
Thomas White's avatar
Thomas White committed
343
344
345
346
347
348
349
	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
350
351
352
353
354
	if ( beam == NULL ) {
		ERROR("You must provide a beam parameters file.\n");
		return 1;
	}

355
	if ( reference_file != NULL ) {
Thomas White's avatar
Thomas White committed
356

357
		RefList *list;
358

359
		list = read_reflections(reference_file);
360
		free(reference_file);
361
		if ( list == NULL ) return 1;
362
		reference_list = asymmetric_indices(list, sym);
363
		reflist_free(list);
364
365
		reference = intensities_from_list(reference_list);

366
367
	}

Thomas White's avatar
Thomas White committed
368
	n_total_patterns = count_patterns(fh);
Thomas White's avatar
Thomas White committed
369
370
371
372
	if ( n_total_patterns == 0 ) {
		ERROR("No patterns to process.\n");
		return 1;
	}
Thomas White's avatar
Thomas White committed
373
374
	STATUS("There are %i patterns to process\n", n_total_patterns);

375
376
	gsl_set_error_handler_off();

377
378
379
380
381
382
383
384
	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);
385
	scalable = new_items();
386
	nobs = 0;
387
388
	for ( i=0; i<n_total_patterns; i++ ) {

389
		RefList *as;
Thomas White's avatar
Thomas White committed
390
		struct image *cur;
391

Thomas White's avatar
Thomas White committed
392
		cur = &images[n_usable_patterns];
393

Thomas White's avatar
Thomas White committed
394
395
396
		cur->det = det;

		if ( read_chunk(fh, cur) != 0 ) {
Thomas White's avatar
Thomas White committed
397
398
			/* Should not happen, because we counted the patterns
			 * earlier. */
Thomas White's avatar
Thomas White committed
399
			ERROR("Failed to read chunk from the input stream.\n");
400
401
402
			return 1;
		}

403
		/* Won't be needing this, if it exists */
Thomas White's avatar
Thomas White committed
404
405
		image_feature_list_free(cur->features);
		cur->features = NULL;
406

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

Thomas White's avatar
Thomas White committed
410
		/* Fill in initial estimates of stuff */
Thomas White's avatar
Thomas White committed
411
412
413
414
415
		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
416
		cur->profile_radius = 0.003e9;
Thomas White's avatar
Thomas White committed
417
		cur->pr_dud = 0;
418

Thomas White's avatar
Thomas White committed
419
		/* Muppet proofing */
Thomas White's avatar
Thomas White committed
420
421
422
		cur->data = NULL;
		cur->flags = NULL;
		cur->beam = NULL;
423

Thomas White's avatar
Thomas White committed
424
		/* This is the raw list of reflections */
Thomas White's avatar
Thomas White committed
425
426
427
		as = asymmetric_indices(cur->reflections, sym);
		reflist_free(cur->reflections);
		cur->reflections = as;
428

429
430
		predict_corresponding_reflections(cur, sym, &n_expected,
		                                 &n_found, &n_notfound);
431

432
433
		nobs += select_scalable_reflections(cur->reflections, scalable);

434
		progress_bar(i, n_total_patterns-1, "Loading pattern data");
Thomas White's avatar
Thomas White committed
435
		n_usable_patterns++;
436
437
438

	}
	fclose(fh);
Thomas White's avatar
Thomas White committed
439
440
	STATUS("Found %5.2f%% of the expected peaks (missed %i of %i).\n",
	       100.0 * (double)n_found / n_expected, n_notfound, n_expected);
441
	STATUS("Mean measurements per scalable unique reflection: %5.2f\n",
442
	       (double)nobs / num_items(scalable));
443

Thomas White's avatar
Thomas White committed
444
	cref = find_common_reflections(images, n_usable_patterns);
Thomas White's avatar
Thomas White committed
445

446
	/* Make initial estimates */
447
	STATUS("Performing initial scaling.\n");
448
	full = scale_intensities(images, n_usable_patterns, sym,
449
	                         scalable, cref, reference);
450

451
452
	select_reflections_for_refinement(images, n_usable_patterns, scalable);

Thomas White's avatar
Thomas White committed
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
	for ( i=0; i<num_items(scalable); i++ ) {
		Reflection *f;
		struct refl_item *it = get_item(scalable, i);
		f = find_refl(full, it->h, it->k, it->l);
		if ( f == NULL ) {
			ERROR("%3i %3i %3i was designated scalable, but no"
			      " full intensity was recorded.\n",
			      it->h, it->k, it->l);
		}
	}

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

		Reflection *refl;
		RefListIterator *iter;

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

			if ( !get_scalable(refl) ) continue;
			get_indices(refl, &h, &k, &l);

			if ( find_item(scalable, h, k, l) == 0 ) {
				ERROR("%3i %3i %3i in image %i is scalable"
				      " but is not in the list of scalable"
				      " reflections.\n", h, k, l, i);
			}
		}
Thomas White's avatar
Thomas White committed
484

Thomas White's avatar
Thomas White committed
485
486
	}

Thomas White's avatar
Thomas White committed
487
	/* Iterate */
488
	for ( i=0; i<n_iter; i++ ) {
Thomas White's avatar
Thomas White committed
489

490
491
		FILE *fhg;
		FILE *fhp;
492
		char filename[1024];
493
		int j;
494

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

497
498
499
500
501
502
503
504
505
506
		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 ) {
507
508
509
510
			ERROR("Failed to open '%s'\n", filename);
			/* Nothing will be written later */
		}

511
512
		if ( reference == NULL ) reference_list = full;

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

517
518
		nobs = 0;
		clear_items(scalable);
519
520
521
522
523
524
525
		for ( j=0; j<n_usable_patterns; j++ ) {

			struct image *cur = &images[j];

			predict_corresponding_reflections(cur, sym, &n_expected,
		                                          &n_found,
		                                          &n_notfound);
526
527
528
529
530
531
532
533

			nobs += select_scalable_reflections(cur->reflections,
			                                    scalable);

		}
		STATUS("Mean measurements per scalable unique "
		       "reflection: %5.2f\n", (double)nobs/num_items(scalable));

Thomas White's avatar
Thomas White committed
534
		/* Re-estimate all the full intensities */
535
		reflist_free(full);
Thomas White's avatar
Thomas White committed
536
		full = scale_intensities(images, n_usable_patterns,
537
		                         sym, scalable, cref, reference);
Thomas White's avatar
Thomas White committed
538

539
540
		fclose(fhg);
		fclose(fhp);
Thomas White's avatar
Thomas White committed
541
542
543
544

	}

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

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

Thomas White's avatar
Thomas White committed
555
	/* Clean up */
Thomas White's avatar
Thomas White committed
556
	for ( i=0; i<n_usable_patterns; i++ ) {
Thomas White's avatar
Thomas White committed
557
		reflist_free(images[i].reflections);
Thomas White's avatar
Thomas White committed
558
	}
559
	reflist_free(full);
560
	delete_items(scalable);
Thomas White's avatar
Thomas White committed
561
	free(sym);
Thomas White's avatar
Thomas White committed
562
	free(outfile);
Thomas White's avatar
Thomas White committed
563
	free_detector_geometry(det);
564
	free(beam);
565
	free(cref);
566
567
568
569
	if ( reference != NULL ) {
		free(reference);
		reflist_free(reference_list);
	}
Thomas White's avatar
Thomas White committed
570
	for ( i=0; i<n_usable_patterns; i++ ) {
571
572
573
574
		cell_free(images[i].indexed_cell);
		free(images[i].filename);
	}
	free(images);
575
576
577

	return 0;
}