process_hkl.c 12.4 KB
Newer Older
1
2
3
4
5
/*
 * process_hkl.c
 *
 * Assemble and process FEL Bragg intensities
 *
Thomas White's avatar
Thomas White committed
6
 * (c) 2006-2010 Thomas White <taw@physics.org>
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
 *
 * 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>

#include "utils.h"
#include "statistics.h"
#include "sfac.h"
27
#include "reflections.h"
28
#include "likelihood.h"
29
#include "symmetry.h"
30
31


Thomas White's avatar
Thomas White committed
32
33
34
35
/* Number of divisions for R vs |q| graphs */
#define RVQDV (20)


36
37
38
static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
Thomas White's avatar
Thomas White committed
39
40
41
	printf(
"Assemble and process FEL Bragg intensities.\n"
"\n"
42
43
44
45
"  -h, --help                Display this help message.\n"
"  -i, --input=<filename>    Specify input filename (\"-\" for stdin).\n"
"  -o, --output=<filename>   Specify output filename for merged intensities\n"
"                             (don't specify for no output).\n"
46
"  -p, --pdb=<filename>      PDB file to use (default: molecule.pdb).\n"
Thomas White's avatar
Thomas White committed
47
"\n"
48
49
50
51
52
53
54
55
56
57
58
59
"      --max-only            Take the integrated intensity to be equal to the\n"
"                             maximum intensity measured for that reflection.\n"
"                             The default is to use the mean value from all\n"
"                             measurements.\n"
"      --sum                 Sum (rather than average) the intensities for the\n"
"                             final output list.  This is useful for comparing\n"
"                             results to radially summed powder patterns, but\n"
"                             will break R-factor analysis.\n"
"      --stop-after=<n>      Stop after processing n patterns.  Zero means\n"
"                             keep going until the end of the input, and is\n"
"                             the default.\n"
"\n"
60
61
"      --scale               Scale each pattern for best fit with the current\n"
"                             model.\n"
Thomas White's avatar
Thomas White committed
62
"  -y, --symmetry=<sym>      Merge according to point group <sym>.\n"
63
);
64
65
66
}


67
68
69
70
71
/* Note "holo" needn't actually be a holohedral point group, if you want to try
 * something strange like resolving from a low-symmetry group into an even
 * lower symmetry one.
 */
static ReflItemList *get_twin_possibilities(const char *holo, const char *mero)
72
{
73
	ReflItemList *test_items;
74
	ReflItemList *twins;
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
	int np;

	np = num_general_equivs(holo) / num_general_equivs(mero);

	test_items = new_items();

	/* Some arbitrarily chosen reflections which can't be special
	 * reflections in any point group, i.e. lots of odd numbers,
	 * prime numbers and so on.  There's probably an analytical
	 * way of working these out, but this will do. */
	add_item(test_items, 1, 2, 3);
	add_item(test_items, 3, 7, 13);
	add_item(test_items, 5, 2, 1);

	twins = get_twins(test_items, holo, mero);
	delete_items(test_items);

Thomas White's avatar
Thomas White committed
92
93
	/* Idiot check.  Wouldn't be necessary if I could prove that the above
	 * set of arbitrarily chosen reflections were always general. */
94
95
96
97
98
99
100
101
102
	if ( num_items(twins) != np ) {
		ERROR("Whoops! Couldn't find all the twinning possiblities.\n");
		abort();
	}

	return twins;
}


Thomas White's avatar
Thomas White committed
103
104
105
static int resolve_twin(const double *model, ReflItemList *observed,
                        const double *patt, ReflItemList *items,
                        ReflItemList *twins, const char *holo, const char *mero)
106
{
107
108
109
110
111
112
113
114
115
116
117
118
119
	int n, i;
	double best_fom = 0.0;
	int best_op = 0;

	n = num_items(twins);

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

		int j;
		int op;
		double *trial_ints = new_list_intensity();
		unsigned int *trial_counts = new_list_count();
		double fom;
Thomas White's avatar
Thomas White committed
120
		ReflItemList *intersection;
121
122
123
124
125
126
127
128
129
130

		op = get_item(twins, i)->op;

		for ( j=0; j<num_items(items); j++ ) {

			signed int h, k, l;
			struct refl_item *r = get_item(items, j);

			get_general_equiv(r->h, r->k, r->l, &h, &k, &l,
			                  holo, op);
131
			get_asymm(h, k, l, &h, &k, &l, mero);
132
133
134
135
136
137
138

			set_intensity(trial_ints, h, k, l,
			              lookup_intensity(patt, r->h, r->k, r->l));
			set_count(trial_counts, h, k, l, 1);

		}

Thomas White's avatar
Thomas White committed
139
140
141
		intersection = intersection_items(observed, items);
		fom = stat_pearson(trial_ints, model, intersection);
		delete_items(intersection);
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158

		free(trial_ints);
		free(trial_counts);

		//printf(" %f", fom);
		if ( fom > best_fom ) {
			best_fom = fom;
			best_op = op;
		}

	}
	//printf("\n");

	return best_op;
}


Thomas White's avatar
Thomas White committed
159
160
161
static void merge_pattern(double *model, ReflItemList *observed,
                          const double *new,  ReflItemList *items,
                          unsigned int *model_counts,  int mo,
162
163
                          ReflItemList *twins,
                          const char *holo, const char *mero)
164
{
165
	int i;
166
	int twin;
167
	ReflItemList *sym_items = new_items();
168

169
	if ( twins != NULL ) {
Thomas White's avatar
Thomas White committed
170
171
		twin = resolve_twin(model, observed, new, items,
		                    twins, holo, mero);
172
173
174
	} else {
		twin = 0;
	}
175

176
	for ( i=0; i<num_items(items); i++ ) {
177
178

		double intensity;
179
		signed int hs, ks, ls;
180
181
		signed int h, k, l;
		struct refl_item *item;
182

183
184
		item = get_item(items, i);

185
186
187
188
		hs = item->h;
		ks = item->k;
		ls = item->l;

189
190
		/* Transform into correct side of the twin law.
		 * "twin" is always zero if no de-twinning is performed. */
191
		get_general_equiv(hs, ks, ls, &h, &k, &l, holo, twin);
192
193

		/* Put into the asymmetric cell for the target group */
194
		get_asymm(h, k, l, &h, &k, &l, mero);
195
196
197

		intensity = lookup_intensity(new, h, k, l);

Thomas White's avatar
Thomas White committed
198
		/* User asked for max only? */
199
200
201
202
203
204
205
206
		if ( !mo ) {
			integrate_intensity(model, h, k, l, intensity);
		} else {
			if ( intensity > lookup_intensity(model, h, k, l) ) {
				set_intensity(model, h, k, l, intensity);
			}
		}

207
208
209
210
211
212
213
214
215
216
		/* Already seen this reflection in this pattern? Complain. */
		if ( !find_item(sym_items, h, k, l) ) {
			/* Add the asymmetric version of this reflection to our
			 * temporary list.  One reflection (in the asymmetric
			 * unit) may appear more than once per pattern if
			 * symmetrically related reflections are present.
			 * That's fine... */
		}	add_item(sym_items, h, k, l);

		/* Increase count count */
Thomas White's avatar
Thomas White committed
217
218
		integrate_count(model_counts, h, k, l, 1);

219
	}
220
221
222
223
224

	/* Dump the reflections in this pattern into the overall list */
	union_items(observed, sym_items);

	delete_items(sym_items);
225
226
227
}


Thomas White's avatar
Thomas White committed
228
229
static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved,
                      unsigned int **pcounts,
230
231
232
233
234
235
236
237
238
239
240
241
                      int config_maxonly, int config_scale, int config_sum,
                      int config_stopafter,
                      ReflItemList *twins, const char *holo, const char *mero,
                      int n_total_patterns)
{
	char *rval;
	float f0;
	int n_nof0 = 0;
	int f0_valid = 0;
	int n_patterns = 0;
	double *new_pattern = new_list_intensity();
	ReflItemList *items = new_items();
Thomas White's avatar
Thomas White committed
242
243
244
245
	ReflItemList *observed = new_items();
	double *model = new_list_intensity();
	unsigned int *counts = new_list_count();
	int i;
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271

	do {

		char line[1024];
		signed int h, k, l;
		float intensity;
		int r;

		rval = fgets(line, 1023, fh);
		if ( (strncmp(line, "Reflections from indexing", 25) == 0)
		    || (strncmp(line, "New pattern", 11) == 0) ) {

			/* Start of first pattern? */
			if ( n_patterns == 0 ) {
				n_patterns++;
				continue;
			}

			/* Assume a default I0 if we don't have one by now */
			if ( config_scale && !f0_valid ) {
				n_nof0++;
				f0 = 1.0;
			}

			/* Scale if requested */
			if ( config_scale ) {
Thomas White's avatar
Thomas White committed
272
273
274
				scale_intensities(model, observed,
				                  new_pattern, items,
				                  f0, f0_valid);
275
276
277
			}

			/* Start of second or later pattern */
Thomas White's avatar
Thomas White committed
278
279
280
			merge_pattern(model, observed, new_pattern, items,
			              counts, config_maxonly,
			              twins, holo, mero);
281
282
283

			if ( n_patterns == config_stopafter ) break;

284
			/* Reset for the next pattern */
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
			n_patterns++;
			clear_items(items);

			progress_bar(n_patterns, n_total_patterns, "Merging");

			f0_valid = 0;

		}

		if ( strncmp(line, "f0 = ", 5) == 0 ) {
			r = sscanf(line, "f0 = %f", &f0);
			if ( r != 1 ) {
				f0 = 1.0;
				f0_valid = 0;
				continue;
			}
			f0_valid = 1;
		}

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

307
		/* Not interested in the central beam */
308
309
		if ( (h==0) && (k==0) && (l==0) ) continue;

310
311
		/* The same raw indices (before mapping into the asymmetric
		 * unit should not turn up twice in one pattern. */
312
313
314
315
316
		if ( find_item(items, h, k, l) != 0 ) {
			ERROR("More than one measurement for %i %i %i in"
			      " pattern number %i\n", h, k, l, n_patterns);
		}
		set_intensity(new_pattern, h, k, l, intensity);
317
318
319

		/* NB: This list contains raw indices, before working out
		 * where they belong in the asymmetric unit. */
320
321
322
323
324
325
326
		add_item(items, h, k, l);

	} while ( rval != NULL );

	delete_items(items);
	free(new_pattern);

327
328
	/* Calculate mean intensity if necessary */
	if ( !config_sum && !config_maxonly ) {
Thomas White's avatar
Thomas White committed
329
330
331
332
333
334
335
336
337
338
339
		for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
			if ( counts[i] > 0 ) {
				model[i] /= (double)counts[i];
			}
		}
	}

	*pmodel = model;
	*pcounts = counts;
	*pobserved = observed;

340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
	STATUS("%i patterns had no f0 valid value.\n", n_nof0);
}


static int count_patterns(FILE *fh)
{
	char *rval;

	int n_total_patterns = 0;
	do {
		char line[1024];

		rval = fgets(line, 1023, fh);
		if ( (strncmp(line, "Reflections from indexing", 25) == 0)
		    || (strncmp(line, "New pattern", 11) == 0) ) {
		    n_total_patterns++;
		}
	} while ( rval != NULL );

	return n_total_patterns;
}


363
364
365
366
int main(int argc, char *argv[])
{
	int c;
	char *filename = NULL;
367
	char *output = NULL;
368
	FILE *fh;
369
	double *model;
Thomas White's avatar
Thomas White committed
370
	unsigned int *counts;
371
	UnitCell *cell;
Thomas White's avatar
Thomas White committed
372
	int config_maxonly = 0;
Thomas White's avatar
Thomas White committed
373
	int config_stopafter = 0;
374
	int config_sum = 0;
375
	int config_scale = 0;
376
	unsigned int n_total_patterns;
377
	char *sym = NULL;
378
	char *pdb = NULL;
379
	ReflItemList *twins;
Thomas White's avatar
Thomas White committed
380
	ReflItemList *observed;
381
	int i;
382
	const char *holo = NULL;
383
384
385
386
387

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
388
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
389
		{"max-only",           0, &config_maxonly,     1},
Thomas White's avatar
Thomas White committed
390
		{"output-every",       1, NULL,               'e'},
Thomas White's avatar
Thomas White committed
391
		{"stop-after",         1, NULL,               's'},
392
		{"sum",                0, &config_sum,         1},
393
		{"scale",              0, &config_scale,       1},
394
		{"symmetry",           1, NULL,               'y'},
395
		{"pdb",                1, NULL,               'p'},
396
397
398
399
		{0, 0, NULL, 0}
	};

	/* Short options */
400
401
	while ((c = getopt_long(argc, argv, "hi:e:ro:p:y:",
	                        longopts, NULL)) != -1) {
402
403

		switch (c) {
Thomas White's avatar
Thomas White committed
404
		case 'h' :
405
406
407
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
408
		case 'i' :
409
410
411
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
412
		case 'o' :
413
414
415
			output = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
416
		case 's' :
Thomas White's avatar
Thomas White committed
417
418
419
			config_stopafter = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
420
		case 'p' :
421
422
423
			pdb = strdup(optarg);
			break;

424
425
426
427
		case 'y' :
			sym = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
428
		case 0 :
429
430
			break;

Thomas White's avatar
Thomas White committed
431
		default :
432
433
434
435
436
437
438
439
440
441
			return 1;
		}

	}

	if ( filename == NULL ) {
		ERROR("Please specify filename using the -i option\n");
		return 1;
	}

442
443
444
445
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

Thomas White's avatar
Thomas White committed
446
447
448
	cell = load_cell_from_pdb(pdb);
	free(pdb);

449
	/* Show useful symmetry information */
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
	if ( sym != NULL ) {
		holo = get_holohedral(sym);
		int np = num_general_equivs(holo) / num_general_equivs(sym);
		if ( np > 1 ) {

			STATUS("Resolving point group %s into %s "
			       "(%i possibilities)\n",
			       holo, sym, np);
			/* Get the list of twin/Bijvoet possibilities */
			twins = get_twin_possibilities(holo, sym);
			STATUS("Twin/inversion operation indices from %s are:",
			       holo);
			for ( i=0; i<num_items(twins); i++ ) {
				STATUS(" %i", get_item(twins, i)->op);
			}
			STATUS("\n");
466

467
468
469
470
		} else {
			STATUS("No twin/inversion resolution necessary.\n");
			twins = NULL;
		}
471
	} else {
472
		STATUS("Not performing any twin/inversion resolution.\n");
473
		twins = NULL;
474
475
		sym = strdup("1");
		holo = strdup("1");
476
	}
477

478
479
480
481
482
483
484
485
486
487
488
	/* Open the data stream */
	if ( strcmp(filename, "-") == 0 ) {
		fh = stdin;
	} else {
		fh = fopen(filename, "r");
	}
	free(filename);
	if ( fh == NULL ) {
		ERROR("Failed to open input file\n");
		return 1;
	}
Thomas White's avatar
Tidy-up    
Thomas White committed
489

490
491
492
493
	/* Count the number of patterns in the file */
	n_total_patterns = count_patterns(fh);
	STATUS("There are %i patterns to process\n", n_total_patterns);
	rewind(fh);
494

Thomas White's avatar
Thomas White committed
495
	merge_all(fh, &model, &observed, &counts,
496
497
498
	          config_maxonly, config_scale, config_sum, config_stopafter,
                  twins, holo, sym, n_total_patterns);
	rewind(fh);
499
500
501

	fclose(fh);

502
	if ( output != NULL ) {
Thomas White's avatar
Thomas White committed
503
		write_reflections(output, observed, model, NULL, counts, cell);
504
	}
Thomas White's avatar
Thomas White committed
505

506
	free(sym);
Thomas White's avatar
Thomas White committed
507
	free(model);
Thomas White's avatar
Thomas White committed
508
	free(counts);
Thomas White's avatar
Thomas White committed
509
510
	free(output);
	free(cell);
511

512
513
	return 0;
}