partialator.c 10.3 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
30
#include "reflections.h"
#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"
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
67
68
	const char *sym;
	ReflItemList *obs;
	double *i_full;
	struct image *image;
69
	FILE *graph;
70
	FILE *pgraph;
Thomas White's avatar
Thomas White committed
71
72
73
74
75
76
77
};


static void refine_image(int mytask, void *tasks)
{
	struct refine_args *all_args = tasks;
	struct refine_args *pargs = &all_args[mytask];
78
79
	struct image *image = pargs->image;

80
	pr_refine(image, pargs->i_full, pargs->sym);
81
82
83
}


84
85
86
87
static void refine_all(struct image *images, int n_total_patterns,
                       struct detector *det, const char *sym,
                       ReflItemList *obs, double *i_full, int nthreads,
                       FILE *graph, FILE *pgraph)
88
{
89
90
	struct refine_args *tasks;
	int i;
Thomas White's avatar
Thomas White committed
91

92
93
94
95
96
97
98
99
100
101
102
	tasks = malloc(n_total_patterns * sizeof(struct refine_args));
	for ( i=0; i<n_total_patterns; i++ ) {

		tasks[i].sym = sym;
		tasks[i].obs = obs;
		tasks[i].i_full = i_full;
		tasks[i].image = &images[i];
		tasks[i].graph = graph;
		tasks[i].pgraph = pgraph;

	}
Thomas White's avatar
Thomas White committed
103

104
105
106
107
108
109
110
	run_thread_range(n_total_patterns, nthreads, "Refining",
	                 refine_image, tasks);

	free(tasks);
}


111
112
113
114
/* Decide which reflections can be scaled */
static void select_scalable_reflections(struct image *images, int n)
{
	int m;
115
	int n_scalable = 0;
116
117
118

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

Thomas White's avatar
Thomas White committed
119
		Reflection *refl;
Thomas White's avatar
Thomas White committed
120
		RefListIterator *iter;
121

Thomas White's avatar
Thomas White committed
122
		for ( refl = first_refl(images[m].reflections, &iter);
Thomas White's avatar
Thomas White committed
123
		      refl != NULL;
Thomas White's avatar
Thomas White committed
124
		      refl = next_refl(refl, iter) ) {
125
126

			int scalable = 1;
Thomas White's avatar
Thomas White committed
127
			double v;
128

Thomas White's avatar
Thomas White committed
129
130
131
			if ( get_partiality(refl) < 0.1 ) scalable = 0;
			v = fabs(get_intensity(refl));
			if ( v < 0.1 ) scalable = 0;
132

Thomas White's avatar
Thomas White committed
133
			set_scalable(refl, scalable);
134
			if ( scalable ) n_scalable++;
135
136
137
138

		}

	}
139
	STATUS("%i reflections selected as scalable.\n", n_scalable);
140
141
142
}


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

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

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

		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;

200
201
202
203
		case 'y' :
			sym = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
204
205
206
207
		case 'o' :
			outfile = strdup(optarg);
			break;

208
209
210
211
		case 'n' :
			n_iter = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
212
213
214
215
216
217
218
219
220
		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
221
222
223
224
225
226
227
228
229
		case 0 :
			break;

		default :
			return 1;
		}

	}

Thomas White's avatar
Thomas White committed
230
	/* Sanitise input filename and open */
Thomas White's avatar
Thomas White committed
231
232
233
234
235
236
237
238
239
240
241
242
243
244
	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
245
246
247
248
249
	/* Sanitise output filename */
	if ( outfile == NULL ) {
		outfile = strdup("facetron.hkl");
	}

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

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

Thomas White's avatar
Thomas White committed
265
266
267
	n_total_patterns = count_patterns(fh);
	STATUS("There are %i patterns to process\n", n_total_patterns);

268
269
270
271
272
273
274
275
	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
276
	obs = new_items();
277
278
279
280
	for ( i=0; i<n_total_patterns; i++ ) {

		UnitCell *cell;
		char *filename;
281
282
283
284
285
		char *rval;
		char line[1024];
		RefList *peaks;
		RefList *transfer;
		Reflection *refl;
Thomas White's avatar
Thomas White committed
286
		RefListIterator *iter;
287
		double ph_en;
288

289
290
		if ( find_chunk(fh, &cell, &filename, &ph_en) == 1 ) {
			ERROR("Couldn't get all of the filenames, cells and"
Thomas White's avatar
Thomas White committed
291
			      " wavelengths from the input stream.\n");
292
293
294
295
			return 1;
		}

		images[i].indexed_cell = cell;
296
		images[i].filename = filename;
Thomas White's avatar
Thomas White committed
297
298
		images[i].div = beam->divergence;
		images[i].bw = beam->bandwidth;
299
		images[i].det = det;
300
		images[i].osf = 1.0;
301
		images[i].profile_radius = 0.005e9;
302
		images[i].reflections = reflist_new();
303
		images[i].lambda = ph_en_to_lambda(eV_to_J(ph_en));
304

Thomas White's avatar
Thomas White committed
305
306
307
		/* Muppet proofing */
		images[i].data = NULL;
		images[i].flags = NULL;
308
		images[i].beam = NULL;
Thomas White's avatar
Thomas White committed
309

310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
		/* Read integrated intensities from pattern */
		peaks = reflist_new();
		do {

			Reflection *refl;
			signed int h, k, l;
			float intensity;
			int r;

			rval = fgets(line, 1023, fh);
			chomp(line);

			if ( (strlen(line) == 0) || (rval == NULL) ) break;

			r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity);
			if ( r != 4 ) continue;

			refl = add_refl(peaks, h, k, l);
			set_int(refl, intensity);

		} while ( (strlen(line) != 0) && (rval != NULL) );

		/* Calculate initial partialities and fill in intensities from
		 * the stream */
		transfer = find_intersections(&images[i], cell, 0);
		images[i].reflections = reflist_new();

Thomas White's avatar
Thomas White committed
337
		for ( refl = first_refl(transfer, &iter);
338
		      refl != NULL;
Thomas White's avatar
Thomas White committed
339
		      refl = next_refl(refl, iter) ) {
340
341
342
343
344
345
346
347

			Reflection *peak;
			Reflection *new;
			signed int h, k, l, ha, ka, la;
			double r1, r2, p, x, y;
			int clamp1, clamp2;

			get_indices(refl, &h, &k, &l);
348
349
			get_detector_pos(refl, &x, &y);
			n_expected++;
Thomas White's avatar
Thomas White committed
350

351
352
			peak = find_refl(peaks, h, k, l);
			if ( peak == NULL ) {
353
				n_notfound++;
354
355
				continue;
			}
356
			n_found++;
357
358

			get_asymm(h, k, l, &ha, &ka, &la, sym);
Thomas White's avatar
Thomas White committed
359
360
361
			if ( find_item(obs, ha, ka, la) == 0 ) {
				add_item(obs, ha, ka, la);
			}
362
363
364
			new = add_refl(images[i].reflections, ha, ka, la);
			get_partial(refl, &r1, &r2, &p, &clamp1, &clamp2);
			get_detector_pos(refl, &x, &y);
365
			set_int(new, get_intensity(peak));
366
367
368
369
370
371
			set_partial(new, r1, r2, p, clamp1, clamp2);
			set_detector_pos(new, 0.0, x, y);

		}
		reflist_free(peaks);
		reflist_free(transfer);
372
		optimise_reflist(images[i].reflections);
373

374
375
376
377
		progress_bar(i, n_total_patterns-1, "Loading pattern data");

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

383
	cref = find_common_reflections(images, n_total_patterns);
Thomas White's avatar
Thomas White committed
384
385
	cts = new_list_count();

386
	/* Make initial estimates */
387
	STATUS("Performing initial scaling.\n");
388
	select_scalable_reflections(images, n_total_patterns);
389
	I_full = scale_intensities(images, n_total_patterns, sym, obs, cref);
390

Thomas White's avatar
Thomas White committed
391
	/* Iterate */
392
	for ( i=0; i<n_iter; i++ ) {
Thomas White's avatar
Thomas White committed
393

394
395
		FILE *fhg;
		FILE *fhp;
396
397
		char filename[1024];

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

400
401
402
403
404
405
406
407
408
409
		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 ) {
410
411
412
413
			ERROR("Failed to open '%s'\n", filename);
			/* Nothing will be written later */
		}

Thomas White's avatar
Thomas White committed
414
		/* Refine the geometry of all patterns to get the best fit */
Thomas White's avatar
Thomas White committed
415
		refine_all(images, n_total_patterns, det, sym, obs, I_full,
416
		           nthreads, fhg, fhp);
Thomas White's avatar
Thomas White committed
417
418

		/* Re-estimate all the full intensities */
Thomas White's avatar
Thomas White committed
419
		free(I_full);
420
		select_scalable_reflections(images, n_total_patterns);
421
422
		I_full = scale_intensities(images, n_total_patterns,
		                           sym, obs, cref);
Thomas White's avatar
Thomas White committed
423

424
425
		fclose(fhg);
		fclose(fhp);
Thomas White's avatar
Thomas White committed
426
427
428
429
430

	}

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

	/* Output results */
Thomas White's avatar
Thomas White committed
435
	write_reflections(outfile, obs, I_full, NULL, NULL, cts, NULL);
Thomas White's avatar
Thomas White committed
436

Thomas White's avatar
Thomas White committed
437
	/* Clean up */
Thomas White's avatar
Thomas White committed
438
	for ( i=0; i<n_total_patterns; i++ ) {
Thomas White's avatar
Thomas White committed
439
		reflist_free(images[i].reflections);
Thomas White's avatar
Thomas White committed
440
441
	}
	free(I_full);
Thomas White's avatar
Thomas White committed
442
	delete_items(obs);
Thomas White's avatar
Thomas White committed
443
	free(sym);
Thomas White's avatar
Thomas White committed
444
	free(outfile);
445
446
	free(det->panels);
	free(det);
447
	free(beam);
Thomas White's avatar
Thomas White committed
448
	free(cts);
449
	free(cref);
450
451
452
453
454
	for ( i=0; i<n_total_patterns; i++ ) {
		cell_free(images[i].indexed_cell);
		free(images[i].filename);
	}
	free(images);
455
456
457

	return 0;
}