partialator.c 10.2 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
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
	for ( i=0; i<n_total_patterns; i++ ) {

279
280
		RefList *predicted;
		RefList *measured;
281
		Reflection *refl;
Thomas White's avatar
Thomas White committed
282
		RefListIterator *iter;
283

284
		if ( read_chunk(fh, &images[i]) == 1 ) {
Thomas White's avatar
Thomas White committed
285
			ERROR("Failed to read chunk from the input stream.\n");
286
287
288
			return 1;
		}

289
290
291
292
		/* 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
293
294
		images[i].div = beam->divergence;
		images[i].bw = beam->bandwidth;
295
		images[i].det = det;
296
297
		images[i].width = det->max_fs;
		images[i].height = det->max_ss;
298
		images[i].osf = 1.0;
299
		images[i].profile_radius = 0.005e9;
300

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

306
307
		/* Calculate initial partialities and fill in intensities from
		 * the stream */
308
309
310
311
312
313
		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;
314
315
		images[i].reflections = reflist_new();

316
		for ( refl = first_refl(predicted, &iter);
317
		      refl != NULL;
Thomas White's avatar
Thomas White committed
318
		      refl = next_refl(refl, iter) ) {
319

320
			Reflection *peak_in_pattern;
321
322
323
324
325
			Reflection *new;
			signed int h, k, l, ha, ka, la;
			double r1, r2, p, x, y;
			int clamp1, clamp2;

326
			/* Get predicted indices and location */
327
			get_indices(refl, &h, &k, &l);
328
329
			get_detector_pos(refl, &x, &y);
			n_expected++;
Thomas White's avatar
Thomas White committed
330

331
332
333
			/* Look for this reflection in the pattern */
			peak_in_pattern = find_refl(measured, h, k, l);
			if ( peak_in_pattern == NULL ) {
334
				n_notfound++;
335
336
				continue;
			}
337
			n_found++;
338

339
			/* Put it into the asymmetric cell */
340
			get_asymm(h, k, l, &ha, &ka, &la, sym);
Thomas White's avatar
Thomas White committed
341
342
343
			if ( find_item(obs, ha, ka, la) == 0 ) {
				add_item(obs, ha, ka, la);
			}
344
345

			/* Create new reflection and copy data across */
346
347
348
			new = add_refl(images[i].reflections, ha, ka, la);
			get_partial(refl, &r1, &r2, &p, &clamp1, &clamp2);
			get_detector_pos(refl, &x, &y);
349
			set_int(new, get_intensity(peak_in_pattern));
350
351
352
353
			set_partial(new, r1, r2, p, clamp1, clamp2);
			set_detector_pos(new, 0.0, x, y);

		}
354
355
356
357
		reflist_free(measured);
		reflist_free(predicted);

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

360
361
362
363
		progress_bar(i, n_total_patterns-1, "Loading pattern data");

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

369
	cref = find_common_reflections(images, n_total_patterns);
Thomas White's avatar
Thomas White committed
370
371
	cts = new_list_count();

372
	/* Make initial estimates */
373
	STATUS("Performing initial scaling.\n");
374
	select_scalable_reflections(images, n_total_patterns);
375
	I_full = scale_intensities(images, n_total_patterns, sym, obs, cref);
376

Thomas White's avatar
Thomas White committed
377
	/* Iterate */
378
	for ( i=0; i<n_iter; i++ ) {
Thomas White's avatar
Thomas White committed
379

380
381
		FILE *fhg;
		FILE *fhp;
382
383
		char filename[1024];

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

386
387
388
389
390
391
392
393
394
395
		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 ) {
396
397
398
399
			ERROR("Failed to open '%s'\n", filename);
			/* Nothing will be written later */
		}

Thomas White's avatar
Thomas White committed
400
		/* Refine the geometry of all patterns to get the best fit */
Thomas White's avatar
Thomas White committed
401
		refine_all(images, n_total_patterns, det, sym, obs, I_full,
402
		           nthreads, fhg, fhp);
Thomas White's avatar
Thomas White committed
403
404

		/* Re-estimate all the full intensities */
Thomas White's avatar
Thomas White committed
405
		free(I_full);
406
		select_scalable_reflections(images, n_total_patterns);
407
408
		I_full = scale_intensities(images, n_total_patterns,
		                           sym, obs, cref);
Thomas White's avatar
Thomas White committed
409

410
411
		fclose(fhg);
		fclose(fhp);
Thomas White's avatar
Thomas White committed
412
413
414
415
416

	}

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

	/* Output results */
Thomas White's avatar
Thomas White committed
421
	//write_reflist(outfile, obs, I_full, NULL, NULL, cts, NULL);
Thomas White's avatar
Thomas White committed
422

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

	return 0;
}