partialator.c 12.1 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;
}


Thomas White's avatar
Thomas White committed
194
195
196
197
int main(int argc, char *argv[])
{
	int c;
	char *infile = NULL;
Thomas White's avatar
Thomas White committed
198
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
199
	char *geomfile = NULL;
Thomas White's avatar
Thomas White committed
200
201
	char *sym = NULL;
	FILE *fh;
Thomas White's avatar
Thomas White committed
202
203
	int nthreads = 1;
	struct detector *det;
204
	ReflItemList *scalable;
Thomas White's avatar
Thomas White committed
205
206
	int i;
	int n_total_patterns;
207
208
	struct image *images;
	int n_iter = 10;
Thomas White's avatar
Thomas White committed
209
	struct beam_params *beam = NULL;
210
	RefList *full;
211
212
213
	int n_found = 0;
	int n_expected = 0;
	int n_notfound = 0;
214
	char *cref;
Thomas White's avatar
Thomas White committed
215
	int n_usable_patterns = 0;
216
	int nobs;
217
	char *reference_file = NULL;
218
	double *reference = NULL;
219
	RefList *reference_list = NULL;
220
	int n_dud;
Thomas White's avatar
Thomas White committed
221
222
223
224
225

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
226
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
227
		{"geometry",           1, NULL,               'g'},
Thomas White's avatar
Thomas White committed
228
		{"beam",               1, NULL,               'b'},
229
		{"symmetry",           1, NULL,               'y'},
230
		{"iterations",         1, NULL,               'n'},
231
		{"reference",          1, NULL,                1},
Thomas White's avatar
Thomas White committed
232
233
234
235
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
236
	while ((c = getopt_long(argc, argv, "hi:g:x:j:y:o:b:",
237
238
	                        longopts, NULL)) != -1)
	{
Thomas White's avatar
Thomas White committed
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256

		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;

257
258
259
260
		case 'y' :
			sym = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
261
262
263
264
		case 'o' :
			outfile = strdup(optarg);
			break;

265
266
267
268
		case 'n' :
			n_iter = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
269
270
271
272
273
274
275
276
277
		case 'b' :
			beam = get_beam_parameters(optarg);
			if ( beam == NULL ) {
				ERROR("Failed to load beam parameters"
				      " from '%s'\n", optarg);
				return 1;
			}
			break;

278
		case 1 :
279
			reference_file = strdup(optarg);
280
281
			break;

Thomas White's avatar
Thomas White committed
282
283
284
285
286
287
288
289
290
		case 0 :
			break;

		default :
			return 1;
		}

	}

Thomas White's avatar
Thomas White committed
291
	/* Sanitise input filename and open */
Thomas White's avatar
Thomas White committed
292
293
294
295
296
297
298
299
300
301
302
303
304
305
	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
306
307
	/* Sanitise output filename */
	if ( outfile == NULL ) {
Thomas White's avatar
Whoops    
Thomas White committed
308
		outfile = strdup("partialator.hkl");
Thomas White's avatar
Thomas White committed
309
310
	}

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

Thomas White's avatar
Thomas White committed
313
	/* Get detector geometry */
Thomas White's avatar
Thomas White committed
314
315
316
317
318
319
320
	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
321
322
323
324
325
	if ( beam == NULL ) {
		ERROR("You must provide a beam parameters file.\n");
		return 1;
	}

326
	if ( reference_file != NULL ) {
Thomas White's avatar
Thomas White committed
327

328
		RefList *list;
329

330
		list = read_reflections(reference_file);
331
		free(reference_file);
332
		if ( list == NULL ) return 1;
333
		reference_list = asymmetric_indices(list, sym);
334
		reflist_free(list);
335
336
		reference = intensities_from_list(reference_list);

337
338
	}

Thomas White's avatar
Thomas White committed
339
	n_total_patterns = count_patterns(fh);
Thomas White's avatar
Thomas White committed
340
341
342
343
	if ( n_total_patterns == 0 ) {
		ERROR("No patterns to process.\n");
		return 1;
	}
Thomas White's avatar
Thomas White committed
344
345
	STATUS("There are %i patterns to process\n", n_total_patterns);

346
347
	gsl_set_error_handler_off();

348
349
350
351
352
353
354
355
	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);
356
	scalable = new_items();
357
	nobs = 0;
358
359
	for ( i=0; i<n_total_patterns; i++ ) {

360
		RefList *as;
Thomas White's avatar
Thomas White committed
361
		struct image *cur;
362

Thomas White's avatar
Thomas White committed
363
		cur = &images[n_usable_patterns];
364

Thomas White's avatar
Thomas White committed
365
366
367
		cur->det = det;

		if ( read_chunk(fh, cur) != 0 ) {
Thomas White's avatar
Thomas White committed
368
369
			/* Should not happen, because we counted the patterns
			 * earlier. */
Thomas White's avatar
Thomas White committed
370
			ERROR("Failed to read chunk from the input stream.\n");
371
372
373
			return 1;
		}

374
		/* Won't be needing this, if it exists */
Thomas White's avatar
Thomas White committed
375
376
		image_feature_list_free(cur->features);
		cur->features = NULL;
377

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

Thomas White's avatar
Thomas White committed
381
		/* Fill in initial estimates of stuff */
Thomas White's avatar
Thomas White committed
382
383
384
385
386
387
388
		cur->div = beam->divergence;
		cur->bw = beam->bandwidth;
		cur->width = det->max_fs;
		cur->height = det->max_ss;
		cur->osf = 1.0;
		cur->profile_radius = 0.0001e9;
		cur->pr_dud = 0;
389

Thomas White's avatar
Thomas White committed
390
		/* Muppet proofing */
Thomas White's avatar
Thomas White committed
391
392
393
		cur->data = NULL;
		cur->flags = NULL;
		cur->beam = NULL;
394

Thomas White's avatar
Thomas White committed
395
		/* This is the raw list of reflections */
Thomas White's avatar
Thomas White committed
396
397
398
		as = asymmetric_indices(cur->reflections, sym);
		reflist_free(cur->reflections);
		cur->reflections = as;
399

400
		update_partialities(cur, sym,
401
		                    &n_expected, &n_found, &n_notfound);
402

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

405
		progress_bar(i, n_total_patterns-1, "Loading pattern data");
Thomas White's avatar
Thomas White committed
406
		n_usable_patterns++;
407
408
409

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

Thomas White's avatar
Thomas White committed
415
	cref = find_common_reflections(images, n_usable_patterns);
Thomas White's avatar
Thomas White committed
416

417
	/* Make initial estimates */
418
	STATUS("Performing initial scaling.\n");
419
	full = scale_intensities(images, n_usable_patterns, sym,
420
	                         scalable, cref, reference);
421

Thomas White's avatar
Thomas White committed
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
	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
453

Thomas White's avatar
Thomas White committed
454
455
	}

Thomas White's avatar
Thomas White committed
456
	/* Iterate */
457
	for ( i=0; i<n_iter; i++ ) {
Thomas White's avatar
Thomas White committed
458

459
460
		FILE *fhg;
		FILE *fhp;
461
462
		char filename[1024];

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

465
466
467
468
469
470
471
472
473
474
		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 ) {
475
476
477
478
			ERROR("Failed to open '%s'\n", filename);
			/* Nothing will be written later */
		}

479
480
		if ( reference == NULL ) reference_list = full;

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

485
486
487
488
489
490
491
492
493
494
495
496
		nobs = 0;
		clear_items(scalable);
		for ( i=0; i<n_usable_patterns; i++ ) {

			struct image *cur = &images[i];
			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
497
		/* Re-estimate all the full intensities */
498
		reflist_free(full);
Thomas White's avatar
Thomas White committed
499
		full = scale_intensities(images, n_usable_patterns,
500
		                         sym, scalable, cref, reference);
Thomas White's avatar
Thomas White committed
501

502
503
		fclose(fhg);
		fclose(fhp);
Thomas White's avatar
Thomas White committed
504
505
506
507

	}

	STATUS("Final scale factors:\n");
508
	n_dud = 0;
Thomas White's avatar
Thomas White committed
509
	for ( i=0; i<n_usable_patterns; i++ ) {
510
		if ( images[i].pr_dud ) n_dud++;
Thomas White's avatar
Thomas White committed
511
		STATUS("%4i : %5.2f\n", i, images[i].osf);
Thomas White's avatar
Thomas White committed
512
	}
513
	STATUS("%i images could not be refined on the last cycle.\n", n_dud);
Thomas White's avatar
Thomas White committed
514
515

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

Thomas White's avatar
Thomas White committed
518
	/* Clean up */
Thomas White's avatar
Thomas White committed
519
	for ( i=0; i<n_usable_patterns; i++ ) {
Thomas White's avatar
Thomas White committed
520
		reflist_free(images[i].reflections);
Thomas White's avatar
Thomas White committed
521
	}
522
	reflist_free(full);
523
	delete_items(scalable);
Thomas White's avatar
Thomas White committed
524
	free(sym);
Thomas White's avatar
Thomas White committed
525
	free(outfile);
Thomas White's avatar
Thomas White committed
526
	free_detector_geometry(det);
527
	free(beam);
528
	free(cref);
529
530
531
532
	if ( reference != NULL ) {
		free(reference);
		reflist_free(reference_list);
	}
Thomas White's avatar
Thomas White committed
533
	for ( i=0; i<n_usable_patterns; i++ ) {
534
535
536
537
		cell_free(images[i].indexed_cell);
		free(images[i].filename);
	}
	free(images);
538
539
540

	return 0;
}