partialator.c 10.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++;

120
	progress_bar(qargs->n_done+1, 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
146

	run_threads(nthreads, refine_image, get_image, done_image,
147
	            &qargs, n_total_patterns-1, 0, 0, 0);
148
149
150
}


Thomas White's avatar
Thomas White committed
151
152
153
154
int main(int argc, char *argv[])
{
	int c;
	char *infile = NULL;
Thomas White's avatar
Thomas White committed
155
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
156
	char *geomfile = NULL;
Thomas White's avatar
Thomas White committed
157
158
	char *sym = NULL;
	FILE *fh;
Thomas White's avatar
Thomas White committed
159
160
	int nthreads = 1;
	struct detector *det;
161
	ReflItemList *scalable;
Thomas White's avatar
Thomas White committed
162
163
	int i;
	int n_total_patterns;
164
165
	struct image *images;
	int n_iter = 10;
Thomas White's avatar
Thomas White committed
166
	struct beam_params *beam = NULL;
167
	RefList *full;
168
169
170
	int n_found = 0;
	int n_expected = 0;
	int n_notfound = 0;
171
	char *cref;
Thomas White's avatar
Thomas White committed
172
	int n_usable_patterns = 0;
173
174
	char *reference_file = NULL;
	RefList *reference;
Thomas White's avatar
Thomas White committed
175
176
177
178
179

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
180
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
181
		{"geometry",           1, NULL,               'g'},
Thomas White's avatar
Thomas White committed
182
		{"beam",               1, NULL,               'b'},
183
		{"symmetry",           1, NULL,               'y'},
184
		{"iterations",         1, NULL,               'n'},
185
		{"reference",          1, NULL,                1},
Thomas White's avatar
Thomas White committed
186
187
188
189
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
190
	while ((c = getopt_long(argc, argv, "hi:g:x:j:y:o:b:",
191
192
	                        longopts, NULL)) != -1)
	{
Thomas White's avatar
Thomas White committed
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210

		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;

211
212
213
214
		case 'y' :
			sym = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
215
216
217
218
		case 'o' :
			outfile = strdup(optarg);
			break;

219
220
221
222
		case 'n' :
			n_iter = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
223
224
225
226
227
228
229
230
231
		case 'b' :
			beam = get_beam_parameters(optarg);
			if ( beam == NULL ) {
				ERROR("Failed to load beam parameters"
				      " from '%s'\n", optarg);
				return 1;
			}
			break;

232
233
234
235
		case 1 :
			reference = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
236
237
238
239
240
241
242
243
244
		case 0 :
			break;

		default :
			return 1;
		}

	}

Thomas White's avatar
Thomas White committed
245
	/* Sanitise input filename and open */
Thomas White's avatar
Thomas White committed
246
247
248
249
250
251
252
253
254
255
256
257
258
259
	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
260
261
	/* Sanitise output filename */
	if ( outfile == NULL ) {
Thomas White's avatar
Whoops    
Thomas White committed
262
		outfile = strdup("partialator.hkl");
Thomas White's avatar
Thomas White committed
263
264
	}

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

Thomas White's avatar
Thomas White committed
267
	/* Get detector geometry */
Thomas White's avatar
Thomas White committed
268
269
270
271
272
273
274
	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
275
276
277
278
279
	if ( beam == NULL ) {
		ERROR("You must provide a beam parameters file.\n");
		return 1;
	}

280
281
282
283
284
285
	if ( reference_file != NULL ) {
		reference = read_reflections(reference_file);
		free(reference_file);
		if ( reference == NULL ) return 1;
	}

Thomas White's avatar
Thomas White committed
286
	n_total_patterns = count_patterns(fh);
Thomas White's avatar
Thomas White committed
287
288
289
290
	if ( n_total_patterns == 0 ) {
		ERROR("No patterns to process.\n");
		return 1;
	}
Thomas White's avatar
Thomas White committed
291
292
	STATUS("There are %i patterns to process\n", n_total_patterns);

293
294
	gsl_set_error_handler_off();

295
296
297
298
299
300
301
302
	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);
303
	scalable = new_items();
304
305
	for ( i=0; i<n_total_patterns; i++ ) {

306
307
		RefList *as;

Thomas White's avatar
Thomas White committed
308
		images[n_usable_patterns].det = NULL;
309

Thomas White's avatar
Thomas White committed
310
		if ( read_chunk(fh, &images[n_usable_patterns]) != 0 ) {
Thomas White's avatar
Thomas White committed
311
312
			/* Should not happen, because we counted the patterns
			 * earlier. */
Thomas White's avatar
Thomas White committed
313
			ERROR("Failed to read chunk from the input stream.\n");
314
315
316
			return 1;
		}

317
		/* Won't be needing this, if it exists */
Thomas White's avatar
Thomas White committed
318
319
		image_feature_list_free(images[n_usable_patterns].features);
		images[n_usable_patterns].features = NULL;
320

Thomas White's avatar
Thomas White committed
321
322
		/* "n_usable_patterns" will not be incremented in this case */
		if ( images[n_usable_patterns].indexed_cell == NULL ) continue;
323

Thomas White's avatar
Thomas White committed
324
325
326
327
328
329
330
331
		/* Fill in initial estimates of stuff */
		images[n_usable_patterns].div = beam->divergence;
		images[n_usable_patterns].bw = beam->bandwidth;
		images[n_usable_patterns].det = det;
		images[n_usable_patterns].width = det->max_fs;
		images[n_usable_patterns].height = det->max_ss;
		images[n_usable_patterns].osf = 1.0;
		images[n_usable_patterns].profile_radius = 0.005e9;
332

Thomas White's avatar
Thomas White committed
333
334
335
336
		/* Muppet proofing */
		images[n_usable_patterns].data = NULL;
		images[n_usable_patterns].flags = NULL;
		images[n_usable_patterns].beam = NULL;
337

Thomas White's avatar
Thomas White committed
338
		/* This is the raw list of reflections */
339
340
341
342
343
		as = asymmetric_indices(images[n_usable_patterns].reflections,
		                        sym);
		optimise_reflist(as);
		reflist_free(images[n_usable_patterns].reflections);
		images[n_usable_patterns].reflections = as;
344

345
346
		update_partialities(&images[n_usable_patterns], sym, scalable,
		                    &n_expected, &n_found, &n_notfound);
347

348
		progress_bar(i, n_total_patterns-1, "Loading pattern data");
Thomas White's avatar
Thomas White committed
349
		n_usable_patterns++;
350
351
352

	}
	fclose(fh);
Thomas White's avatar
Thomas White committed
353
354
	STATUS("Found %5.2f%% of the expected peaks (missed %i of %i).\n",
	       100.0 * (double)n_found / n_expected, n_notfound, n_expected);
355
356
	STATUS("Mean measurements per scalable unique reflection: %5.2f\n",
	       (double)n_found / num_items(scalable));
357

Thomas White's avatar
Thomas White committed
358
	cref = find_common_reflections(images, n_usable_patterns);
Thomas White's avatar
Thomas White committed
359

360
	/* Make initial estimates */
361
	STATUS("Performing initial scaling.\n");
362
363
	full = scale_intensities(images, n_usable_patterns, sym,
	                         scalable, cref);
364

Thomas White's avatar
Thomas White committed
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
	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
398
	/* Iterate */
399
	for ( i=0; i<n_iter; i++ ) {
Thomas White's avatar
Thomas White committed
400

401
402
		FILE *fhg;
		FILE *fhp;
403
404
		char filename[1024];

405
		STATUS("Post refinement iteration %i of %i\n", i+1, n_iter);
Thomas White's avatar
Thomas White committed
406

407
408
409
410
411
412
413
414
415
416
		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 ) {
417
418
419
420
			ERROR("Failed to open '%s'\n", filename);
			/* Nothing will be written later */
		}

Thomas White's avatar
Thomas White committed
421
		/* Refine the geometry of all patterns to get the best fit */
422
		refine_all(images, n_total_patterns, det, sym, scalable, full,
423
		           nthreads, fhg, fhp);
Thomas White's avatar
Thomas White committed
424
425

		/* Re-estimate all the full intensities */
426
		reflist_free(full);
Thomas White's avatar
Thomas White committed
427
		full = scale_intensities(images, n_usable_patterns,
428
		                         sym, scalable, cref);
Thomas White's avatar
Thomas White committed
429

430
431
		fclose(fhg);
		fclose(fhp);
Thomas White's avatar
Thomas White committed
432
433
434
435

	}

	STATUS("Final scale factors:\n");
Thomas White's avatar
Thomas White committed
436
	for ( i=0; i<n_usable_patterns; i++ ) {
Thomas White's avatar
Thomas White committed
437
		STATUS("%4i : %5.2f\n", i, images[i].osf);
Thomas White's avatar
Thomas White committed
438
439
440
	}

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

Thomas White's avatar
Thomas White committed
443
	/* Clean up */
Thomas White's avatar
Thomas White committed
444
	for ( i=0; i<n_usable_patterns; i++ ) {
Thomas White's avatar
Thomas White committed
445
		reflist_free(images[i].reflections);
Thomas White's avatar
Thomas White committed
446
	}
447
	reflist_free(full);
448
	delete_items(scalable);
Thomas White's avatar
Thomas White committed
449
	free(sym);
Thomas White's avatar
Thomas White committed
450
	free(outfile);
Thomas White's avatar
Thomas White committed
451
	free_detector_geometry(det);
452
	free(beam);
453
	free(cref);
Thomas White's avatar
Thomas White committed
454
	for ( i=0; i<n_usable_patterns; i++ ) {
455
456
457
458
		cell_free(images[i].indexed_cell);
		free(images[i].filename);
	}
	free(images);
459
460
461

	return 0;
}