partialator.c 9.93 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
115
116
117
/* Decide which reflections can be scaled */
static void select_scalable_reflections(struct image *images, int n)
{
	int m;

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

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

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

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

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

Thomas White's avatar
Thomas White committed
132
			set_scalable(refl, scalable);
133
134
135
136
137
138
139

		}

	}
}


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

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

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

		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;

196
197
198
199
		case 'y' :
			sym = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
200
201
202
203
		case 'o' :
			outfile = strdup(optarg);
			break;

204
205
206
207
		case 'n' :
			n_iter = atoi(optarg);
			break;

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

		default :
			return 1;
		}

	}

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

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

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

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

264
265
266
267
268
269
270
271
	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
272
	obs = new_items();
273
274
275
276
	for ( i=0; i<n_total_patterns; i++ ) {

		UnitCell *cell;
		char *filename;
277
278
279
280
281
		char *rval;
		char line[1024];
		RefList *peaks;
		RefList *transfer;
		Reflection *refl;
Thomas White's avatar
Thomas White committed
282
		RefListIterator *iter;
283
		double ph_en;
284

285
286
		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
287
			      " wavelengths from the input stream.\n");
288
289
290
291
292
			return 1;
		}

		images[i].indexed_cell = cell;

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

Thomas White's avatar
Thomas White committed
302
303
304
		/* Muppet proofing */
		images[i].data = NULL;
		images[i].flags = NULL;
305
		images[i].beam = NULL;
Thomas White's avatar
Thomas White committed
306

307
308
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
		/* 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
334
		for ( refl = first_refl(transfer, &iter);
335
		      refl != NULL;
Thomas White's avatar
Thomas White committed
336
		      refl = next_refl(refl, iter) ) {
337
338
339
340
341
342
343
344

			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);
345
346
			get_detector_pos(refl, &x, &y);
			n_expected++;
Thomas White's avatar
Thomas White committed
347

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

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

		}
		reflist_free(peaks);
		reflist_free(transfer);
368

369
370
371
372
		progress_bar(i, n_total_patterns-1, "Loading pattern data");

	}
	fclose(fh);
373
374
375
	STATUS("Found %5.2f%% of the expected peaks (%i %i %i).\n",
	       100.0 * (double)n_found / n_expected,
	       n_found, n_notfound, n_expected);
376

Thomas White's avatar
Thomas White committed
377
378
	cts = new_list_count();

379
	/* Make initial estimates */
380
	STATUS("Performing initial scaling.\n");
381
	select_scalable_reflections(images, n_total_patterns);
Thomas White's avatar
Thomas White committed
382
	I_full = scale_intensities(images, n_total_patterns, sym, obs);
383

Thomas White's avatar
Thomas White committed
384
	/* Iterate */
385
	for ( i=0; i<n_iter; i++ ) {
Thomas White's avatar
Thomas White committed
386

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

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

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

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

		/* Re-estimate all the full intensities */
Thomas White's avatar
Thomas White committed
412
		free(I_full);
413
		select_scalable_reflections(images, n_total_patterns);
Thomas White's avatar
Thomas White committed
414
		I_full = scale_intensities(images, n_total_patterns, sym, obs);
Thomas White's avatar
Thomas White committed
415

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

	}

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

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

Thomas White's avatar
Thomas White committed
429
	/* Clean up */
Thomas White's avatar
Thomas White committed
430
	for ( i=0; i<n_total_patterns; i++ ) {
Thomas White's avatar
Thomas White committed
431
		reflist_free(images[i].reflections);
Thomas White's avatar
Thomas White committed
432
433
	}
	free(I_full);
Thomas White's avatar
Thomas White committed
434
	delete_items(obs);
Thomas White's avatar
Thomas White committed
435
	free(sym);
Thomas White's avatar
Thomas White committed
436
	free(outfile);
437
438
	free(det->panels);
	free(det);
439
	free(beam);
Thomas White's avatar
Thomas White committed
440
	free(cts);
441
442
443
444
445
	for ( i=0; i<n_total_patterns; i++ ) {
		cell_free(images[i].indexed_cell);
		free(images[i].filename);
	}
	free(images);
446
447
448

	return 0;
}