partialator.c 9.52 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
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;
139
140
141
	/* 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
142
143

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


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

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

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

		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;

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

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

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

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

		default :
			return 1;
		}

	}

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

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

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

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

277
278
279
280
281
282
283
284
	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
285
	obs = new_items();
286
287
	for ( i=0; i<n_total_patterns; i++ ) {

Thomas White's avatar
Thomas White committed
288
		images[n_usable_patterns].det = NULL;
289

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

297
		/* Won't be needing this, if it exists */
Thomas White's avatar
Thomas White committed
298
299
		image_feature_list_free(images[n_usable_patterns].features);
		images[n_usable_patterns].features = NULL;
300

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

Thomas White's avatar
Thomas White committed
304
305
306
307
308
309
310
311
		/* 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;
312

Thomas White's avatar
Thomas White committed
313
314
315
316
		/* Muppet proofing */
		images[n_usable_patterns].data = NULL;
		images[n_usable_patterns].flags = NULL;
		images[n_usable_patterns].beam = NULL;
317

Thomas White's avatar
Thomas White committed
318
319
320
321
		/* This is the raw list of reflections */
		images[n_usable_patterns].raw_reflections =
		                          images[n_usable_patterns].reflections;
		images[n_usable_patterns].reflections = NULL;
322

Thomas White's avatar
Thomas White committed
323
324
325
		update_partialities_and_asymm(&images[n_usable_patterns], sym,
		                              obs, &n_expected, &n_found,
		                              &n_notfound);
326

327
		progress_bar(i, n_total_patterns-1, "Loading pattern data");
Thomas White's avatar
Thomas White committed
328
		n_usable_patterns++;
329
330
331

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

Thomas White's avatar
Thomas White committed
337
	cref = find_common_reflections(images, n_usable_patterns);
Thomas White's avatar
Thomas White committed
338

339
	/* Make initial estimates */
340
	STATUS("Performing initial scaling.\n");
Thomas White's avatar
Thomas White committed
341
	full = scale_intensities(images, n_usable_patterns, sym, obs, cref);
342

Thomas White's avatar
Thomas White committed
343
	/* Iterate */
344
	for ( i=0; i<n_iter; i++ ) {
Thomas White's avatar
Thomas White committed
345

346
347
		FILE *fhg;
		FILE *fhp;
348
349
		char filename[1024];

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

352
353
354
355
356
357
358
359
360
361
		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 ) {
362
363
364
365
			ERROR("Failed to open '%s'\n", filename);
			/* Nothing will be written later */
		}

Thomas White's avatar
Thomas White committed
366
		/* Refine the geometry of all patterns to get the best fit */
367
		refine_all(images, n_total_patterns, det, sym, obs, full,
368
		           nthreads, fhg, fhp);
Thomas White's avatar
Thomas White committed
369
370

		/* Re-estimate all the full intensities */
371
		reflist_free(full);
Thomas White's avatar
Thomas White committed
372
		full = scale_intensities(images, n_usable_patterns,
373
		                         sym, obs, cref);
Thomas White's avatar
Thomas White committed
374

375
376
		fclose(fhg);
		fclose(fhp);
Thomas White's avatar
Thomas White committed
377
378
379
380

	}

	STATUS("Final scale factors:\n");
Thomas White's avatar
Thomas White committed
381
	for ( i=0; i<n_usable_patterns; i++ ) {
Thomas White's avatar
Thomas White committed
382
		STATUS("%4i : %5.2f\n", i, images[i].osf);
Thomas White's avatar
Thomas White committed
383
384
385
	}

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

Thomas White's avatar
Thomas White committed
388
	/* Clean up */
Thomas White's avatar
Thomas White committed
389
	for ( i=0; i<n_usable_patterns; i++ ) {
Thomas White's avatar
Thomas White committed
390
		reflist_free(images[i].reflections);
Thomas White's avatar
Thomas White committed
391
		reflist_free(images[i].raw_reflections);
Thomas White's avatar
Thomas White committed
392
	}
393
	reflist_free(full);
Thomas White's avatar
Thomas White committed
394
	delete_items(obs);
Thomas White's avatar
Thomas White committed
395
	free(sym);
Thomas White's avatar
Thomas White committed
396
	free(outfile);
Thomas White's avatar
Thomas White committed
397
	free_detector_geometry(det);
398
	free(beam);
399
	free(cref);
Thomas White's avatar
Thomas White committed
400
	for ( i=0; i<n_usable_patterns; i++ ) {
401
402
403
404
		cell_free(images[i].indexed_cell);
		free(images[i].filename);
	}
	free(images);
405
406
407

	return 0;
}