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


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

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

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

		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;

216
217
218
219
		case 'y' :
			sym = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
220
221
222
223
		case 'o' :
			outfile = strdup(optarg);
			break;

224
225
226
227
		case 'n' :
			n_iter = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
228
229
230
231
232
233
234
235
236
		case 'b' :
			beam = get_beam_parameters(optarg);
			if ( beam == NULL ) {
				ERROR("Failed to load beam parameters"
				      " from '%s'\n", optarg);
				return 1;
			}
			break;

237
		case 1 :
238
			reference_file = strdup(optarg);
239
240
			break;

Thomas White's avatar
Thomas White committed
241
242
243
244
245
246
247
248
249
		case 0 :
			break;

		default :
			return 1;
		}

	}

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

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

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

285
	if ( reference_file != NULL ) {
Thomas White's avatar
Thomas White committed
286

287
		RefList *list;
288

289
		list = read_reflections(reference_file);
290
		free(reference_file);
291
		if ( list == NULL ) return 1;
292
		reference_list = asymmetric_indices(list, sym);
293
		reflist_free(list);
294
295
		reference = intensities_from_list(reference_list);

296
297
	}

Thomas White's avatar
Thomas White committed
298
	n_total_patterns = count_patterns(fh);
Thomas White's avatar
Thomas White committed
299
300
301
302
	if ( n_total_patterns == 0 ) {
		ERROR("No patterns to process.\n");
		return 1;
	}
Thomas White's avatar
Thomas White committed
303
304
	STATUS("There are %i patterns to process\n", n_total_patterns);

305
306
	gsl_set_error_handler_off();

307
308
309
310
311
312
313
314
	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);
315
	scalable = new_items();
316
317
	for ( i=0; i<n_total_patterns; i++ ) {

318
319
		RefList *as;

320
		images[n_usable_patterns].det = det;
321

Thomas White's avatar
Thomas White committed
322
		if ( read_chunk(fh, &images[n_usable_patterns]) != 0 ) {
Thomas White's avatar
Thomas White committed
323
324
			/* Should not happen, because we counted the patterns
			 * earlier. */
Thomas White's avatar
Thomas White committed
325
			ERROR("Failed to read chunk from the input stream.\n");
326
327
328
			return 1;
		}

329
		/* Won't be needing this, if it exists */
Thomas White's avatar
Thomas White committed
330
331
		image_feature_list_free(images[n_usable_patterns].features);
		images[n_usable_patterns].features = NULL;
332

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

Thomas White's avatar
Thomas White committed
336
337
338
339
340
341
		/* Fill in initial estimates of stuff */
		images[n_usable_patterns].div = beam->divergence;
		images[n_usable_patterns].bw = beam->bandwidth;
		images[n_usable_patterns].width = det->max_fs;
		images[n_usable_patterns].height = det->max_ss;
		images[n_usable_patterns].osf = 1.0;
342
		images[n_usable_patterns].profile_radius = 0.0001e9;
343
		images[n_usable_patterns].pr_dud = 0;
344

Thomas White's avatar
Thomas White committed
345
346
347
348
		/* Muppet proofing */
		images[n_usable_patterns].data = NULL;
		images[n_usable_patterns].flags = NULL;
		images[n_usable_patterns].beam = NULL;
349

Thomas White's avatar
Thomas White committed
350
		/* This is the raw list of reflections */
351
352
353
354
355
		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;
356

357
358
		update_partialities(&images[n_usable_patterns], sym, scalable,
		                    &n_expected, &n_found, &n_notfound);
359

360
		progress_bar(i, n_total_patterns-1, "Loading pattern data");
Thomas White's avatar
Thomas White committed
361
		n_usable_patterns++;
362
363
364

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

Thomas White's avatar
Thomas White committed
370
	cref = find_common_reflections(images, n_usable_patterns);
Thomas White's avatar
Thomas White committed
371

372
	/* Make initial estimates */
373
	STATUS("Performing initial scaling.\n");
374
	full = scale_intensities(images, n_usable_patterns, sym,
375
	                         scalable, cref, reference);
376

Thomas White's avatar
Thomas White committed
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
	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
408

Thomas White's avatar
Thomas White committed
409
410
	}

Thomas White's avatar
Thomas White committed
411
	/* Iterate */
412
	for ( i=0; i<n_iter; i++ ) {
Thomas White's avatar
Thomas White committed
413

414
415
		FILE *fhg;
		FILE *fhp;
416
417
		char filename[1024];

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

420
421
422
423
424
425
426
427
428
429
		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 ) {
430
431
432
433
			ERROR("Failed to open '%s'\n", filename);
			/* Nothing will be written later */
		}

434
435
		if ( reference == NULL ) reference_list = full;

Thomas White's avatar
Thomas White committed
436
		/* Refine the geometry of all patterns to get the best fit */
437
438
		refine_all(images, n_usable_patterns, det, sym, scalable,
		           reference_list, nthreads, fhg, fhp);
Thomas White's avatar
Thomas White committed
439
440

		/* Re-estimate all the full intensities */
441
		reflist_free(full);
Thomas White's avatar
Thomas White committed
442
		full = scale_intensities(images, n_usable_patterns,
443
		                         sym, scalable, cref, reference);
Thomas White's avatar
Thomas White committed
444

445
446
		fclose(fhg);
		fclose(fhp);
Thomas White's avatar
Thomas White committed
447
448
449
450

	}

	STATUS("Final scale factors:\n");
451
	n_dud = 0;
Thomas White's avatar
Thomas White committed
452
	for ( i=0; i<n_usable_patterns; i++ ) {
453
		if ( images[i].pr_dud ) n_dud++;
Thomas White's avatar
Thomas White committed
454
		STATUS("%4i : %5.2f\n", i, images[i].osf);
Thomas White's avatar
Thomas White committed
455
	}
456
	STATUS("%i images could not be refined on the last cycle.\n", n_dud);
Thomas White's avatar
Thomas White committed
457
458

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

Thomas White's avatar
Thomas White committed
461
	/* Clean up */
Thomas White's avatar
Thomas White committed
462
	for ( i=0; i<n_usable_patterns; i++ ) {
Thomas White's avatar
Thomas White committed
463
		reflist_free(images[i].reflections);
Thomas White's avatar
Thomas White committed
464
	}
465
	reflist_free(full);
466
	delete_items(scalable);
Thomas White's avatar
Thomas White committed
467
	free(sym);
Thomas White's avatar
Thomas White committed
468
	free(outfile);
Thomas White's avatar
Thomas White committed
469
	free_detector_geometry(det);
470
	free(beam);
471
	free(cref);
472
473
474
475
	if ( reference != NULL ) {
		free(reference);
		reflist_free(reference_list);
	}
Thomas White's avatar
Thomas White committed
476
	for ( i=0; i<n_usable_patterns; i++ ) {
477
478
479
480
		cell_free(images[i].indexed_cell);
		free(images[i].filename);
	}
	free(images);
481
482
483

	return 0;
}