partialator.c 11 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

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


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


Thomas White's avatar
Thomas White committed
63
struct refine_args
Thomas White's avatar
Thomas White committed
64
{
Thomas White's avatar
Thomas White committed
65
66
	const char *sym;
	ReflItemList *obs;
67
	RefList *full;
Thomas White's avatar
Thomas White committed
68
	struct image *image;
69
	FILE *graph;
70
	FILE *pgraph;
Thomas White's avatar
Thomas White committed
71
72
73
};


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

90
	pr_refine(image, pargs->full, pargs->sym);
91
92
93
}


Thomas White's avatar
Thomas White committed
94
static void *get_image(void *vqargs)
95
{
Thomas White's avatar
Thomas White committed
96
97
	struct refine_args *task;
	struct queue_args *qargs = vqargs;
Thomas White's avatar
Thomas White committed
98

Thomas White's avatar
Thomas White committed
99
100
	task = malloc(sizeof(struct refine_args));
	memcpy(task, &qargs->task_defaults, sizeof(struct refine_args));
101

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

Thomas White's avatar
Thomas White committed
104
105
106
107
	qargs->n++;

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

109

Thomas White's avatar
Thomas White committed
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
static void done_image(void *vqargs, void *task)
{
	struct queue_args *qargs = vqargs;

	qargs->n_done++;

	progress_bar(qargs->n_done, qargs->n_total_patterns, "Refining");
	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;
	qargs.n_total_patterns = n_total_patterns;
	qargs.images = images;

	run_threads(nthreads, refine_image, get_image, done_image,
	            &qargs, n_total_patterns, 0, 0, 0);
144
145
146
}


147
148
149
150
/* Decide which reflections can be scaled */
static void select_scalable_reflections(struct image *images, int n)
{
	int m;
151
	int n_scalable = 0;
152
153
154

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

Thomas White's avatar
Thomas White committed
155
		Reflection *refl;
Thomas White's avatar
Thomas White committed
156
		RefListIterator *iter;
157

Thomas White's avatar
Thomas White committed
158
		for ( refl = first_refl(images[m].reflections, &iter);
Thomas White's avatar
Thomas White committed
159
		      refl != NULL;
Thomas White's avatar
Thomas White committed
160
		      refl = next_refl(refl, iter) ) {
161
162

			int scalable = 1;
Thomas White's avatar
Thomas White committed
163
			double v;
164

Thomas White's avatar
Thomas White committed
165
166
167
			if ( get_partiality(refl) < 0.1 ) scalable = 0;
			v = fabs(get_intensity(refl));
			if ( v < 0.1 ) scalable = 0;
168

Thomas White's avatar
Thomas White committed
169
			set_scalable(refl, scalable);
170
			if ( scalable ) n_scalable++;
171
172
173
174

		}

	}
175
	STATUS("%i reflections selected as scalable.\n", n_scalable);
176
177
178
}


Thomas White's avatar
Thomas White committed
179
180
181
182
int main(int argc, char *argv[])
{
	int c;
	char *infile = NULL;
Thomas White's avatar
Thomas White committed
183
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
184
	char *geomfile = NULL;
Thomas White's avatar
Thomas White committed
185
186
	char *sym = NULL;
	FILE *fh;
Thomas White's avatar
Thomas White committed
187
188
	int nthreads = 1;
	struct detector *det;
Thomas White's avatar
Thomas White committed
189
190
191
	ReflItemList *obs;
	int i;
	int n_total_patterns;
192
193
	struct image *images;
	int n_iter = 10;
Thomas White's avatar
Thomas White committed
194
	struct beam_params *beam = NULL;
195
	RefList *full;
196
197
198
	int n_found = 0;
	int n_expected = 0;
	int n_notfound = 0;
199
	char *cref;
Thomas White's avatar
Thomas White committed
200
201
202
203
204

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
205
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
206
		{"geometry",           1, NULL,               'g'},
Thomas White's avatar
Thomas White committed
207
		{"beam",               1, NULL,               'b'},
208
		{"symmetry",           1, NULL,               'y'},
209
		{"iterations",         1, NULL,               'n'},
Thomas White's avatar
Thomas White committed
210
211
212
213
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
214
	while ((c = getopt_long(argc, argv, "hi:g:x:j:y:o:b:",
215
216
	                        longopts, NULL)) != -1)
	{
Thomas White's avatar
Thomas White committed
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234

		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;

235
236
237
238
		case 'y' :
			sym = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
239
240
241
242
		case 'o' :
			outfile = strdup(optarg);
			break;

243
244
245
246
		case 'n' :
			n_iter = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
247
248
249
250
251
252
253
254
255
		case 'b' :
			beam = get_beam_parameters(optarg);
			if ( beam == NULL ) {
				ERROR("Failed to load beam parameters"
				      " from '%s'\n", optarg);
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
256
257
258
259
260
261
262
263
264
		case 0 :
			break;

		default :
			return 1;
		}

	}

Thomas White's avatar
Thomas White committed
265
	/* Sanitise input filename and open */
Thomas White's avatar
Thomas White committed
266
267
268
269
270
271
272
273
274
275
276
277
278
279
	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
280
281
	/* Sanitise output filename */
	if ( outfile == NULL ) {
Thomas White's avatar
Whoops    
Thomas White committed
282
		outfile = strdup("partialator.hkl");
Thomas White's avatar
Thomas White committed
283
284
	}

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

Thomas White's avatar
Thomas White committed
287
	/* Get detector geometry */
Thomas White's avatar
Thomas White committed
288
289
290
291
292
293
294
	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
295
296
297
298
299
	if ( beam == NULL ) {
		ERROR("You must provide a beam parameters file.\n");
		return 1;
	}

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

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);
Thomas White's avatar
Thomas White committed
315
	obs = new_items();
316
317
	for ( i=0; i<n_total_patterns; i++ ) {

318
319
		RefList *predicted;
		RefList *measured;
320
		Reflection *refl;
Thomas White's avatar
Thomas White committed
321
		RefListIterator *iter;
322

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

330
331
332
		/* FIXME: This leaves gaps in the array */
		if ( images[i].indexed_cell == NULL ) continue;

333
334
335
336
		/* Won't be needing this, if it exists */
		image_feature_list_free(images[i].features);
		images[i].features = NULL;

Thomas White's avatar
Thomas White committed
337
338
		images[i].div = beam->divergence;
		images[i].bw = beam->bandwidth;
339
		images[i].det = det;
340
341
		images[i].width = det->max_fs;
		images[i].height = det->max_ss;
342
		images[i].osf = 1.0;
343
		images[i].profile_radius = 0.005e9;
344

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

350
351
		/* Calculate initial partialities and fill in intensities from
		 * the stream */
352
353
354
355
356
357
		predicted = find_intersections(&images[i],
		                               images[i].indexed_cell, 0);

		/* We start again with a new reflection list, this time with
		 * the asymmetric indices */
		measured = images[i].reflections;
358
359
		images[i].reflections = reflist_new();

360
		for ( refl = first_refl(predicted, &iter);
361
		      refl != NULL;
Thomas White's avatar
Thomas White committed
362
		      refl = next_refl(refl, iter) ) {
363

364
			Reflection *peak_in_pattern;
365
366
367
368
369
			Reflection *new;
			signed int h, k, l, ha, ka, la;
			double r1, r2, p, x, y;
			int clamp1, clamp2;

370
			/* Get predicted indices and location */
371
			get_indices(refl, &h, &k, &l);
372
373
			get_detector_pos(refl, &x, &y);
			n_expected++;
Thomas White's avatar
Thomas White committed
374

375
376
377
			/* Look for this reflection in the pattern */
			peak_in_pattern = find_refl(measured, h, k, l);
			if ( peak_in_pattern == NULL ) {
378
				n_notfound++;
379
380
				continue;
			}
381
			n_found++;
382

383
			/* Put it into the asymmetric cell */
384
			get_asymm(h, k, l, &ha, &ka, &la, sym);
Thomas White's avatar
Thomas White committed
385
386
387
			if ( find_item(obs, ha, ka, la) == 0 ) {
				add_item(obs, ha, ka, la);
			}
388
389

			/* Create new reflection and copy data across */
390
391
392
			new = add_refl(images[i].reflections, ha, ka, la);
			get_partial(refl, &r1, &r2, &p, &clamp1, &clamp2);
			get_detector_pos(refl, &x, &y);
393
			set_int(new, get_intensity(peak_in_pattern));
394
395
396
397
			set_partial(new, r1, r2, p, clamp1, clamp2);
			set_detector_pos(new, 0.0, x, y);

		}
398
399
400
401
		reflist_free(measured);
		reflist_free(predicted);

		/* Do magic on the reflection list to make things go faster */
402
		optimise_reflist(images[i].reflections);
403

404
405
406
407
		progress_bar(i, n_total_patterns-1, "Loading pattern data");

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

413
	cref = find_common_reflections(images, n_total_patterns);
Thomas White's avatar
Thomas White committed
414

415
	/* Make initial estimates */
416
	STATUS("Performing initial scaling.\n");
417
	select_scalable_reflections(images, n_total_patterns);
418
	full = scale_intensities(images, n_total_patterns, sym, obs, cref);
419

Thomas White's avatar
Thomas White committed
420
	/* Iterate */
421
	for ( i=0; i<n_iter; i++ ) {
Thomas White's avatar
Thomas White committed
422

423
424
		FILE *fhg;
		FILE *fhp;
425
426
		char filename[1024];

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

429
430
431
432
433
434
435
436
437
438
		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 ) {
439
440
441
442
			ERROR("Failed to open '%s'\n", filename);
			/* Nothing will be written later */
		}

Thomas White's avatar
Thomas White committed
443
		/* Refine the geometry of all patterns to get the best fit */
444
		refine_all(images, n_total_patterns, det, sym, obs, full,
445
		           nthreads, fhg, fhp);
Thomas White's avatar
Thomas White committed
446
447

		/* Re-estimate all the full intensities */
448
		reflist_free(full);
449
		select_scalable_reflections(images, n_total_patterns);
450
451
		full = scale_intensities(images, n_total_patterns,
		                         sym, obs, cref);
Thomas White's avatar
Thomas White committed
452

453
454
		fclose(fhg);
		fclose(fhp);
Thomas White's avatar
Thomas White committed
455
456
457
458
459

	}

	STATUS("Final scale factors:\n");
	for ( i=0; i<n_total_patterns; i++ ) {
Thomas White's avatar
Thomas White committed
460
		STATUS("%4i : %5.2f\n", i, images[i].osf);
Thomas White's avatar
Thomas White committed
461
462
463
	}

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

Thomas White's avatar
Thomas White committed
466
	/* Clean up */
Thomas White's avatar
Thomas White committed
467
	for ( i=0; i<n_total_patterns; i++ ) {
Thomas White's avatar
Thomas White committed
468
		reflist_free(images[i].reflections);
Thomas White's avatar
Thomas White committed
469
	}
470
	reflist_free(full);
Thomas White's avatar
Thomas White committed
471
	delete_items(obs);
Thomas White's avatar
Thomas White committed
472
	free(sym);
Thomas White's avatar
Thomas White committed
473
	free(outfile);
474
475
	free(det->panels);
	free(det);
476
	free(beam);
477
	free(cref);
478
479
480
481
482
	for ( i=0; i<n_total_patterns; i++ ) {
		cell_free(images[i].indexed_cell);
		free(images[i].filename);
	}
	free(images);
483
484
485

	return 0;
}