partialator.c 10.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
"\n"
Thomas White's avatar
Thomas White committed
60
61
62
63
"  -j <n>                     Run <n> analyses in parallel.\n");
}


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


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

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


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

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

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

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

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

110

Thomas White's avatar
Thomas White committed
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
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;
140
141
142
	/* FIXME: Not refining the first image, for now */
	qargs.n_total_patterns = n_total_patterns-1;
	qargs.images = images+1;
Thomas White's avatar
Thomas White committed
143
144

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


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

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

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

		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;

206
207
208
209
		case 'y' :
			sym = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
210
211
212
213
		case 'o' :
			outfile = strdup(optarg);
			break;

214
215
216
217
		case 'n' :
			n_iter = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
218
219
220
221
222
223
224
225
226
		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
227
228
229
230
231
232
233
234
235
		case 0 :
			break;

		default :
			return 1;
		}

	}

Thomas White's avatar
Thomas White committed
236
	/* Sanitise input filename and open */
Thomas White's avatar
Thomas White committed
237
238
239
240
241
242
243
244
245
246
247
248
249
250
	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
251
252
	/* Sanitise output filename */
	if ( outfile == NULL ) {
Thomas White's avatar
Whoops    
Thomas White committed
253
		outfile = strdup("partialator.hkl");
Thomas White's avatar
Thomas White committed
254
255
	}

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

Thomas White's avatar
Thomas White committed
258
	/* Get detector geometry */
Thomas White's avatar
Thomas White committed
259
260
261
262
263
264
265
	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
266
267
268
269
270
	if ( beam == NULL ) {
		ERROR("You must provide a beam parameters file.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
271
	n_total_patterns = count_patterns(fh);
Thomas White's avatar
Thomas White committed
272
273
274
275
	if ( n_total_patterns == 0 ) {
		ERROR("No patterns to process.\n");
		return 1;
	}
Thomas White's avatar
Thomas White committed
276
277
	STATUS("There are %i patterns to process\n", n_total_patterns);

278
279
	gsl_set_error_handler_off();

280
281
282
283
284
285
286
287
	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);
288
	scalable = new_items();
289
290
	for ( i=0; i<n_total_patterns; i++ ) {

291
292
		RefList *as;

Thomas White's avatar
Thomas White committed
293
		images[n_usable_patterns].det = NULL;
294

Thomas White's avatar
Thomas White committed
295
		if ( read_chunk(fh, &images[n_usable_patterns]) != 0 ) {
Thomas White's avatar
Thomas White committed
296
297
			/* Should not happen, because we counted the patterns
			 * earlier. */
Thomas White's avatar
Thomas White committed
298
			ERROR("Failed to read chunk from the input stream.\n");
299
300
301
			return 1;
		}

302
		/* Won't be needing this, if it exists */
Thomas White's avatar
Thomas White committed
303
304
		image_feature_list_free(images[n_usable_patterns].features);
		images[n_usable_patterns].features = NULL;
305

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

Thomas White's avatar
Thomas White committed
309
310
311
312
313
314
315
316
		/* 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;
317

Thomas White's avatar
Thomas White committed
318
319
320
321
		/* Muppet proofing */
		images[n_usable_patterns].data = NULL;
		images[n_usable_patterns].flags = NULL;
		images[n_usable_patterns].beam = NULL;
322

Thomas White's avatar
Thomas White committed
323
		/* This is the raw list of reflections */
324
325
326
327
328
		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;
329

330
331
		update_partialities(&images[n_usable_patterns], sym, scalable,
		                    &n_expected, &n_found, &n_notfound);
332

333
		progress_bar(i, n_total_patterns-1, "Loading pattern data");
Thomas White's avatar
Thomas White committed
334
		n_usable_patterns++;
335
336
337

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

Thomas White's avatar
Thomas White committed
343
	cref = find_common_reflections(images, n_usable_patterns);
Thomas White's avatar
Thomas White committed
344

345
	/* Make initial estimates */
346
	STATUS("Performing initial scaling.\n");
347
348
	full = scale_intensities(images, n_usable_patterns, sym,
	                         scalable, cref);
349

Thomas White's avatar
Thomas White committed
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
	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
383
	/* Iterate */
384
	for ( i=0; i<n_iter; i++ ) {
Thomas White's avatar
Thomas White committed
385

386
387
		FILE *fhg;
		FILE *fhp;
388
389
		char filename[1024];

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

392
393
394
395
396
397
398
399
400
401
		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 ) {
402
403
404
405
			ERROR("Failed to open '%s'\n", filename);
			/* Nothing will be written later */
		}

Thomas White's avatar
Thomas White committed
406
		/* Refine the geometry of all patterns to get the best fit */
407
		refine_all(images, n_total_patterns, det, sym, scalable, full,
408
		           nthreads, fhg, fhp);
Thomas White's avatar
Thomas White committed
409
410

		/* Re-estimate all the full intensities */
411
		reflist_free(full);
Thomas White's avatar
Thomas White committed
412
		full = scale_intensities(images, n_usable_patterns,
413
		                         sym, scalable, cref);
Thomas White's avatar
Thomas White committed
414

415
416
		fclose(fhg);
		fclose(fhp);
Thomas White's avatar
Thomas White committed
417
418
419
420

	}

	STATUS("Final scale factors:\n");
Thomas White's avatar
Thomas White committed
421
	for ( i=0; i<n_usable_patterns; i++ ) {
Thomas White's avatar
Thomas White committed
422
		STATUS("%4i : %5.2f\n", i, images[i].osf);
Thomas White's avatar
Thomas White committed
423
424
425
	}

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

Thomas White's avatar
Thomas White committed
428
	/* Clean up */
Thomas White's avatar
Thomas White committed
429
	for ( i=0; i<n_usable_patterns; i++ ) {
Thomas White's avatar
Thomas White committed
430
		reflist_free(images[i].reflections);
Thomas White's avatar
Thomas White committed
431
	}
432
	reflist_free(full);
433
	delete_items(scalable);
Thomas White's avatar
Thomas White committed
434
	free(sym);
Thomas White's avatar
Thomas White committed
435
	free(outfile);
Thomas White's avatar
Thomas White committed
436
	free_detector_geometry(det);
437
	free(beam);
438
	free(cref);
Thomas White's avatar
Thomas White committed
439
	for ( i=0; i<n_usable_patterns; i++ ) {
440
441
442
443
		cell_free(images[i].indexed_cell);
		free(images[i].filename);
	}
	free(images);
444
445
446

	return 0;
}