get_hkl.c 12.7 KB
Newer Older
1
2
3
/*
 * get_hkl.c
 *
4
 * Small program to manipulate reflection lists
5
 *
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
23
24
 *
 * 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"
Thomas White's avatar
Thomas White committed
25
#include "reflist-utils.h"
26
#include "symmetry.h"
27
#include "beam-parameters.h"
28
29
30
31
32
33


static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
34
"Manipulate reflection lists.\n"
35
"\n"
Thomas White's avatar
Thomas White committed
36
37
"  -h, --help                 Display this help message.\n"
"\n"
38
39
40
41
"  -i, --input=<file>         Read reflections from <file>.\n"
"  -y, --symmetry=<sym>       The symmetry of the input reflection list.\n"
"\n"
"You can add noise to the reflections with either of:\n"
Thomas White's avatar
Thomas White committed
42
"      --poisson              Simulate Poisson samples.\n"
43
"      --noise                Add 10%% random noise.\n"
44
45
46
47
"\n"
"To calculate Poisson samples accurately, you must also give:\n"
"  -b, --beam=<file>          Get beam parameters from file.\n"
"\n"
Thomas White's avatar
Thomas White committed
48
49
"You can artificially 'twin' the reflections, or expand them out.  You can also"
" do both, in which case the 'twinning' will be done first:\n"
50
51
"  -w, --twin=<sym>           Generate twinned data according to the given\n"
"                              point group.\n"
Thomas White's avatar
Thomas White committed
52
"  -e, --expand=<sym>         Expand reflections to this point group.\n"
53
"\n"
54
55
56
57
58
"Use this option with care, and only if you understand why it might sometimes\n"
" be necessary:\n"
"  --trim-centrics            Remove reflections which are duplicated in the\n"
"                              point group specified with the '-y' option.\n"
"\n"
59
60
61
62
"You can restrict which reflections are written out:\n"
"  -t, --template=<filename>  Only include reflections mentioned in file.\n"
"\n"
"You might sometimes need to do this:\n"
63
64
"      --multiplicity         Multiply intensities by the number of\n"
"                              equivalent reflections.\n"
65
66
67
"\n"
"Don't forget to specify the output filename:\n"
"  -o, --output=<filename>    Output filename (default: stdout).\n"
Thomas White's avatar
Thomas White committed
68
69
"  -p, --pdb=<filename>       Use unit cell parameters from this PDB file to\n"
"                              generate resolution values in the output file.\n"
Thomas White's avatar
Thomas White committed
70
);
Thomas White's avatar
Thomas White committed
71
72
73
74
}


/* Apply Poisson noise to all reflections */
Thomas White's avatar
Thomas White committed
75
static void poisson_reflections(RefList *list, double adu_per_photon)
Thomas White's avatar
Thomas White committed
76
{
Thomas White's avatar
Thomas White committed
77
78
	Reflection *refl;
	RefListIterator *iter;
Thomas White's avatar
Thomas White committed
79

Thomas White's avatar
Thomas White committed
80
81
82
	for ( refl = first_refl(list, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
Thomas White's avatar
Thomas White committed
83

Thomas White's avatar
Thomas White committed
84
		double val, c;
Thomas White's avatar
Thomas White committed
85

Thomas White's avatar
Thomas White committed
86
		val = get_intensity(refl);
87

88
		c = adu_per_photon * poisson_noise(val/adu_per_photon);
Thomas White's avatar
Thomas White committed
89
		set_int(refl, c);
Thomas White's avatar
Thomas White committed
90
91

	}
92
93
94
95
}


/* Apply 10% uniform noise to all reflections */
Thomas White's avatar
Thomas White committed
96
static void noise_reflections(RefList *list)
97
{
Thomas White's avatar
Thomas White committed
98
99
	Reflection *refl;
	RefListIterator *iter;
100

Thomas White's avatar
Thomas White committed
101
102
103
	for ( refl = first_refl(list, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
104

Thomas White's avatar
Thomas White committed
105
		double val, r;
106

Thomas White's avatar
Thomas White committed
107
		val = get_intensity(refl);
108
109
110
111

		r = (double)random()/RAND_MAX;
		val += 0.1 * val * r;

Thomas White's avatar
Thomas White committed
112
		set_int(refl, val);
113

Thomas White's avatar
Thomas White committed
114
	}
115
116
117
}


Thomas White's avatar
Thomas White committed
118
static RefList *template_reflections(RefList *list, RefList *template)
119
{
Thomas White's avatar
Thomas White committed
120
121
122
	Reflection *refl;
	RefListIterator *iter;
	RefList *out;
123

Thomas White's avatar
Thomas White committed
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
	out = reflist_new();

	for ( refl = first_refl(template, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {

		signed int h, k, l;
		Reflection *new;
		Reflection *old;

		get_indices(refl, &h, &k, &l);

		old = find_refl(list, h, k, l);
		if ( old == NULL ) continue;

		new = add_refl(out, h, k, l);
		copy_data(new, old);

	}

	return out;
}


static RefList *twin_reflections(RefList *in,
Thomas White's avatar
Thomas White committed
149
                                 const SymOpList *holo, const SymOpList *mero)
Thomas White's avatar
Thomas White committed
150
151
152
153
{
	Reflection *refl;
	RefListIterator *iter;
	RefList *out;
Thomas White's avatar
Thomas White committed
154
	SymOpMask *m;
155

Thomas White's avatar
Thomas White committed
156
	/* FIXME: Check properly by coset decomposition */
Thomas White's avatar
Thomas White committed
157
	if ( num_equivs(holo, NULL) < num_equivs(mero, NULL) ) {
Thomas White's avatar
Thomas White committed
158
159
		ERROR("%s is not a subgroup of %s!\n", symmetry_name(mero),
		                                       symmetry_name(holo));
160
161
162
		return NULL;
	}

Thomas White's avatar
Thomas White committed
163
	out = reflist_new();
Thomas White's avatar
Thomas White committed
164
	m = new_symopmask(holo);
Thomas White's avatar
Thomas White committed
165
166
167
168

	for ( refl = first_refl(in, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
169

170
		double total, sigma;
171
172
173
174
		signed int h, k, l;
		int n, j;
		int skip;

Thomas White's avatar
Thomas White committed
175
		get_indices(refl, &h, &k, &l);
176

177
		/* There is a many-to-one correspondence between reflections
Thomas White's avatar
Thomas White committed
178
179
180
		 * in the merohedral and holohedral groups.  Do the calculation
		 * only once for each reflection in the holohedral group, which
		 * contains fewer reflections.
181
		 */
Thomas White's avatar
Thomas White committed
182
		get_asymm(holo, h, k, l, &h, &k, &l);
Thomas White's avatar
Thomas White committed
183
		if ( find_refl(out, h, k, l) != NULL ) continue;
184

185
186
		total = 0.0;
		sigma = 0.0;
187
		skip = 0;
Thomas White's avatar
Thomas White committed
188
189
190
		special_position(holo, m, h, k, l);
		n = num_equivs(holo, m);

191
192
193
194
195
		for ( j=0; j<n; j++ ) {

			signed int he, ke, le;
			signed int hu, ku, lu;

Thomas White's avatar
Thomas White committed
196
			get_equiv(holo, m, j, h, k, l, &he, &ke, &le);
197
198
199
200
201
202

			/* Do we have this reflection?
			 * We might not have the particular (merohedral)
			 * equivalent which belongs to our definition of the
			 * asymmetric unit cell, so check them all.
			 */
Thomas White's avatar
Thomas White committed
203
204
			if ( !find_equiv_in_list(in, he, ke, le, mero,
			                         &hu, &ku, &lu) ) {
205
206
207
208
209
				/* Don't have this reflection, so bail out */
				ERROR("Twinning %i %i %i requires the %i %i %i "
				      "reflection (or an equivalent in %s), "
				      "which I don't have. %i %i %i won't "
				      "appear in the output\n",
Thomas White's avatar
Thomas White committed
210
211
				      h, k, l, he, ke, le, symmetry_name(mero),
				      h, k, l);
212
213
214
215
				skip = 1;
				break;
			}

Thomas White's avatar
Thomas White committed
216
217
			total += get_intensity(refl);
			sigma += pow(get_esd_intensity(refl), 2.0);
218
219
220
221
222

		}

		if ( !skip ) {

Thomas White's avatar
Thomas White committed
223
224
225
			Reflection *new = add_refl(out, h, k, l);
			set_int(new, total);
			set_esd_intensity(new, sqrt(sigma));
226
227
228
229
230

		}

	}

Thomas White's avatar
Thomas White committed
231
	return out;
232
233
234
}


Thomas White's avatar
Thomas White committed
235
236
static RefList *expand_reflections(RefList *in, const SymOpList *target,
                                                const SymOpList *initial)
Thomas White's avatar
Thomas White committed
237
{
Thomas White's avatar
Thomas White committed
238
239
240
	Reflection *refl;
	RefListIterator *iter;
	RefList *out;
Thomas White's avatar
Thomas White committed
241
	SymOpMask *m;
Thomas White's avatar
Thomas White committed
242

Thomas White's avatar
Thomas White committed
243
	/* FIXME: Check properly */
Thomas White's avatar
Thomas White committed
244
	if ( num_equivs(target, NULL) > num_equivs(initial, NULL) ) {
Thomas White's avatar
Thomas White committed
245
246
		ERROR("%s is not a subgroup of %s!\n", symmetry_name(initial),
		                                       symmetry_name(target));
247
248
249
		return NULL;
	}

Thomas White's avatar
Thomas White committed
250
	out = reflist_new();
Thomas White's avatar
Thomas White committed
251
	m = new_symopmask(initial);
Thomas White's avatar
Thomas White committed
252
253
254
255

	for ( refl = first_refl(in, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
Thomas White's avatar
Thomas White committed
256
257
258
259
260

		signed int h, k, l;
		int n, j;
		double intensity;

Thomas White's avatar
Thomas White committed
261
262
		get_indices(refl, &h, &k, &l);
		intensity = get_intensity(refl);
Thomas White's avatar
Thomas White committed
263

Thomas White's avatar
Thomas White committed
264
265
		special_position(initial, m, h, k, l);
		n = num_equivs(initial, m);
Thomas White's avatar
Thomas White committed
266

Thomas White's avatar
Thomas White committed
267
		/* For each equivalent in the higher symmetry group */
Thomas White's avatar
Thomas White committed
268
269
270
		for ( j=0; j<n; j++ ) {

			signed int he, ke, le;
Thomas White's avatar
Thomas White committed
271
			Reflection *new;
272
273
			int have_phase;
			double ph;
Thomas White's avatar
Thomas White committed
274
275

			/* Get the equivalent */
Thomas White's avatar
Thomas White committed
276
			get_equiv(initial, m, j, h, k, l, &he, &ke, &le);
Thomas White's avatar
Thomas White committed
277
278

			/* Put it into the asymmetric unit for the target */
Thomas White's avatar
Thomas White committed
279
			get_asymm(target, he, ke, le, &he, &ke, &le);
Thomas White's avatar
Thomas White committed
280
281

			/* Make sure the intensity is in the right place */
Thomas White's avatar
Thomas White committed
282
283
			new = add_refl(out, he, ke, le);
			copy_data(new, refl);
Thomas White's avatar
Thomas White committed
284

Thomas White's avatar
Thomas White committed
285
286
			/* FIXME: Make phase negative if the reflection is
			 * separated from the original via an inversion */
287
288
289
290
			get_phase(refl, &have_phase);
			if ( have_phase ) {
				set_ph(new, -ph);
			}
Thomas White's avatar
Thomas White committed
291

Thomas White's avatar
Thomas White committed
292
293
294
295
		}

	}

Thomas White's avatar
Thomas White committed
296
297
	free_symopmask(m);

Thomas White's avatar
Thomas White committed
298
	return out;
Thomas White's avatar
Thomas White committed
299
300
301
}


302
303
304
305
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
334
335
336
337
static RefList *trim_centrics(RefList *in, const SymOpList *sym)
{
	Reflection *refl;
	RefListIterator *iter;
	RefList *out;

	out = reflist_new();

	for ( refl = first_refl(in, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) )
	{
		signed int h, k, l;
		signed int ha, ka, la;
		Reflection *new;

		get_indices(refl, &h, &k, &l);

		/* Put it into the asymmetric unit */
		get_asymm(sym, h, k, l, &ha, &ka, &la);

		new = find_refl(out, ha, ka, la);
		if ( new != NULL ) {
			STATUS("Trimmed %i %i %i\n", h, k, l);
			continue;
		}

		/* Add new reflection under asymmetric (unique) indices */
		new = add_refl(out, ha, ka, la);
		copy_data(new, refl);
	}

	return out;
}


338
339
340
int main(int argc, char *argv[])
{
	int c;
341
342
	int config_noise = 0;
	int config_poisson = 0;
343
	int config_multi = 0;
344
	int config_trimc = 0;
Thomas White's avatar
Thomas White committed
345
346
347
348
349
350
	char *holo_str = NULL;
	char *mero_str = NULL;
	char *expand_str = NULL;
	SymOpList *holo;
	SymOpList *mero;
	SymOpList *expand;
Thomas White's avatar
Thomas White committed
351
352
	char *input_file = NULL;
	char *template = NULL;
Thomas White's avatar
Thomas White committed
353
	char *output = NULL;
354
	char *beamfile = NULL;
355
	struct beam_params *beam = NULL;
Thomas White's avatar
Thomas White committed
356
357
	RefList *input;
	UnitCell *cell = NULL;
358
359
360
361

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
Thomas White's avatar
Thomas White committed
362
		{"template",           1, NULL,               't'},
363
364
		{"poisson",            0, &config_poisson,     1},
		{"noise",              0, &config_noise,       1},
Thomas White's avatar
Thomas White committed
365
		{"output",             1, NULL,               'o'},
366
		{"symmetry",           1, NULL,               'y'},
Thomas White's avatar
Thomas White committed
367
368
		{"twin",               1, NULL,               'w'},
		{"expand",             1, NULL,               'e'},
369
		{"intensities",        1, NULL,               'i'},
370
		{"multiplicity",       0, &config_multi,       1},
371
		{"beam",               1, NULL,               'b'},
Thomas White's avatar
Thomas White committed
372
		{"pdb",                1, NULL,               'p'},
373
		{"trim-centrics",      0, &config_trimc,       1},
374
375
376
377
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
378
	while ((c = getopt_long(argc, argv, "ht:o:i:w:y:e:b:p:",
Thomas White's avatar
Thomas White committed
379
	                        longopts, NULL)) != -1) {
380
381

		switch (c) {
Thomas White's avatar
Thomas White committed
382
		case 'h' :
383
384
385
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
386
		case 't' :
Thomas White's avatar
Thomas White committed
387
388
389
			template = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
390
		case 'o' :
Thomas White's avatar
Thomas White committed
391
392
393
			output = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
394
		case 'i' :
Thomas White's avatar
Thomas White committed
395
			input_file = strdup(optarg);
396
397
			break;

398
		case 'y' :
Thomas White's avatar
Thomas White committed
399
			mero_str = strdup(optarg);
400
401
402
			break;

		case 'w' :
Thomas White's avatar
Thomas White committed
403
			holo_str = strdup(optarg);
404
405
			break;

Thomas White's avatar
Thomas White committed
406
		case 'e' :
Thomas White's avatar
Thomas White committed
407
			expand_str = strdup(optarg);
Thomas White's avatar
Thomas White committed
408
409
			break;

410
411
412
413
		case 'b' :
			beamfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
414
415
416
417
418
419
420
421
		case 'p' :
			cell = load_cell_from_pdb(optarg);
			if ( cell == NULL ) {
				ERROR("Failed to get cell from '%s'\n", optarg);
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
422
		case 0 :
423
424
			break;

Thomas White's avatar
Thomas White committed
425
		default :
426
427
428
429
430
			return 1;
		}

	}

Thomas White's avatar
Thomas White committed
431
	if ( (holo_str != NULL) && (expand_str != NULL) ) {
Thomas White's avatar
Thomas White committed
432
433
		ERROR("You cannot 'twin' and 'expand' at the same time.\n");
		ERROR("Decide which one you want to do first.\n");
434
435
436
437
438
439
440
441
442
443
		return 1;
	}

	if ( beamfile != NULL ) {
		beam = get_beam_parameters(beamfile);
		if ( beam == NULL ) {
			ERROR("Failed to load beam parameters from '%s'\n",
			      beamfile);
			return 1;
		}
Thomas White's avatar
Thomas White committed
444
445
	}

Thomas White's avatar
Thomas White committed
446
447
448
	if ( cell == NULL ) {
		ERROR("You need to give a PDB file with the unit cell.\n");
		return 1;
449
	}
450

Thomas White's avatar
Thomas White committed
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
	if ( holo_str != NULL ) {
		holo = get_pointgroup(holo_str);
		free(holo_str);
	} else {
		holo = NULL;
	}
	if ( mero_str != NULL ) {
		mero = get_pointgroup(mero_str);
		free(mero_str);
	} else {
		mero = NULL;
	}
	if ( expand_str != NULL ) {
		expand = get_pointgroup(expand_str);
		free(expand_str);
	} else {
		expand = NULL;
	}

Thomas White's avatar
Thomas White committed
470
	input = read_reflections(input_file);
Thomas White's avatar
Thomas White committed
471
472
473
474
	if ( input == NULL ) {
		ERROR("Problem reading input file %s\n", input_file);
		return 1;
	}
Thomas White's avatar
Thomas White committed
475
	free(input_file);
476
477
478
479

	STATUS("%i reflections in input.\n", num_reflections(input));

	if ( !config_trimc && check_list_symmetry(input, mero) ) {
480
		ERROR("The input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
481
		      " have symmetry %s\n", symmetry_name(mero));
482
		return 1;
483
	}
484

485
486
	if ( config_poisson ) {
		if ( beam != NULL ) {
Thomas White's avatar
Thomas White committed
487
			poisson_reflections(input, beam->adu_per_photon);
488
489
490
491
492
493
494
		} else {
			ERROR("You must give a beam parameters file in order"
			      " to calculate Poisson noise.\n");
			return 1;
		}
	}

Thomas White's avatar
Thomas White committed
495
	if ( config_noise ) noise_reflections(input);
Thomas White's avatar
Thomas White committed
496

497
	if ( holo != NULL ) {
498

Thomas White's avatar
Thomas White committed
499
		RefList *new;
Thomas White's avatar
Thomas White committed
500
501
		STATUS("Twinning from %s into %s\n", symmetry_name(mero),
		                                     symmetry_name(holo));
Thomas White's avatar
Thomas White committed
502
503
504
505
506
507
508
509
510
		new = twin_reflections(input, holo, mero);

		/* Replace old with new */
		reflist_free(input);
		input = new;

		/* The symmetry of the list has changed */
		free(mero);
		mero = holo;
511

Thomas White's avatar
Thomas White committed
512
513
	}

Thomas White's avatar
Thomas White committed
514
	if ( expand != NULL ) {
515

Thomas White's avatar
Thomas White committed
516
		RefList *new;
Thomas White's avatar
Thomas White committed
517
518
		STATUS("Expanding from %s into %s\n", symmetry_name(mero),
		                                      symmetry_name(expand));
Thomas White's avatar
Thomas White committed
519
520
521
522
523
		new = expand_reflections(input, expand, mero);

		/* Replace old with new */
		reflist_free(input);
		input = new;
524

Thomas White's avatar
Thomas White committed
525
526
	}

527
528
529
530
531
532
533
534
535
536
537
538
539
	if ( config_trimc ) {

		RefList *new;

		STATUS("Trimming duplicate reflections in %s\n",
		       symmetry_name(mero));
		new = trim_centrics(input, mero);
		reflist_free(input);
		input = new;
		STATUS("%i output reflections\n", num_reflections(input));

	}

540
541
	if ( config_multi ) {

Thomas White's avatar
Thomas White committed
542
543
		Reflection *refl;
		RefListIterator *iter;
Thomas White's avatar
Thomas White committed
544
545
546
		SymOpMask *m;

		m = new_symopmask(mero);
547

Thomas White's avatar
Thomas White committed
548
549
550
		for ( refl = first_refl(input, &iter);
		      refl != NULL;
		      refl = next_refl(refl, iter) ) {
551
552

			double inty;
Thomas White's avatar
Thomas White committed
553
554
555
556
			signed int h, k, l;

			get_indices(refl, &h, &k, &l);
			inty = get_intensity(refl);
557

Thomas White's avatar
Thomas White committed
558
559
			special_position(mero, m, h, k, l);
			inty *= (double)num_equivs(mero, m);
Thomas White's avatar
Thomas White committed
560
			set_int(refl, inty);
561
562

		}
Thomas White's avatar
Thomas White committed
563
564

		free_symopmask(m);
565
566
	}

Thomas White's avatar
Thomas White committed
567
	if ( template ) {
568

Thomas White's avatar
Thomas White committed
569
570
571
572
		RefList *t = read_reflections(template);
		RefList *new = template_reflections(input, t);
		reflist_free(input);
		input = new;
573

Thomas White's avatar
Thomas White committed
574
575
	}

Thomas White's avatar
Thomas White committed
576
	write_reflist(output, input, cell);
Thomas White's avatar
Thomas White committed
577

Thomas White's avatar
Thomas White committed
578
	reflist_free(input);
579
580
581

	return 0;
}