get_hkl.c 13.3 KB
Newer Older
1
2
3
/*
 * get_hkl.c
 *
4
 * Small program to manipulate reflection lists
5
 *
Thomas White's avatar
Thomas White committed
6
 * Copyright © 2012 Thomas White <taw@physics.org>
7
 *
Thomas White's avatar
Thomas White committed
8
9
10
11
12
13
14
15
16
17
18
19
20
21
 * This file is part of CrystFEL.
 *
 * CrystFEL is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * CrystFEL is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
 *
 */


#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
38
#include "reflist-utils.h"
39
#include "symmetry.h"
40
#include "beam-parameters.h"
41
42
43
44
45
46


static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
47
"Manipulate reflection lists.\n"
48
"\n"
Thomas White's avatar
Thomas White committed
49
50
"  -h, --help                 Display this help message.\n"
"\n"
51
52
53
54
"  -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
55
"      --poisson              Simulate Poisson samples.\n"
56
"      --noise                Add 10%% random noise.\n"
57
58
59
60
"\n"
"To calculate Poisson samples accurately, you must also give:\n"
"  -b, --beam=<file>          Get beam parameters from file.\n"
"\n"
61
"You can artificially 'twin' the reflections, or expand them out.\n"
62
63
"  -w, --twin=<sym>           Generate twinned data according to the given\n"
"                              point group.\n"
Thomas White's avatar
Thomas White committed
64
"  -e, --expand=<sym>         Expand reflections to this point group.\n"
65
"\n"
66
67
68
69
70
"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"
71
72
73
74
"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"
75
76
"      --multiplicity         Multiply intensities by the number of\n"
"                              equivalent reflections.\n"
77
78
79
"\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
80
81
"  -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
82
);
Thomas White's avatar
Thomas White committed
83
84
85
86
}


/* Apply Poisson noise to all reflections */
Thomas White's avatar
Thomas White committed
87
static void poisson_reflections(RefList *list, double adu_per_photon)
Thomas White's avatar
Thomas White committed
88
{
Thomas White's avatar
Thomas White committed
89
90
	Reflection *refl;
	RefListIterator *iter;
Thomas White's avatar
Thomas White committed
91

Thomas White's avatar
Thomas White committed
92
93
94
	for ( refl = first_refl(list, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
Thomas White's avatar
Thomas White committed
95

Thomas White's avatar
Thomas White committed
96
		double val, c;
Thomas White's avatar
Thomas White committed
97

Thomas White's avatar
Thomas White committed
98
		val = get_intensity(refl);
99

100
		c = adu_per_photon * poisson_noise(val/adu_per_photon);
Thomas White's avatar
Thomas White committed
101
		set_int(refl, c);
Thomas White's avatar
Thomas White committed
102
103

	}
104
105
106
107
}


/* Apply 10% uniform noise to all reflections */
Thomas White's avatar
Thomas White committed
108
static void noise_reflections(RefList *list)
109
{
Thomas White's avatar
Thomas White committed
110
111
	Reflection *refl;
	RefListIterator *iter;
112

Thomas White's avatar
Thomas White committed
113
114
115
	for ( refl = first_refl(list, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
116

Thomas White's avatar
Thomas White committed
117
		double val, r;
118

Thomas White's avatar
Thomas White committed
119
		val = get_intensity(refl);
120
121
122
123

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

Thomas White's avatar
Thomas White committed
124
		set_int(refl, val);
125

Thomas White's avatar
Thomas White committed
126
	}
127
128
129
}


Thomas White's avatar
Thomas White committed
130
static RefList *template_reflections(RefList *list, RefList *template)
131
{
Thomas White's avatar
Thomas White committed
132
133
134
	Reflection *refl;
	RefListIterator *iter;
	RefList *out;
135

Thomas White's avatar
Thomas White committed
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
	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
161
                                 const SymOpList *holo, const SymOpList *mero)
Thomas White's avatar
Thomas White committed
162
163
164
165
{
	Reflection *refl;
	RefListIterator *iter;
	RefList *out;
Thomas White's avatar
Thomas White committed
166
	SymOpMask *m;
167

Thomas White's avatar
Thomas White committed
168
	/* FIXME: Check properly by coset decomposition */
Thomas White's avatar
Thomas White committed
169
	if ( num_equivs(holo, NULL) < num_equivs(mero, NULL) ) {
Thomas White's avatar
Thomas White committed
170
171
		ERROR("%s is not a subgroup of %s!\n", symmetry_name(mero),
		                                       symmetry_name(holo));
172
173
174
		return NULL;
	}

Thomas White's avatar
Thomas White committed
175
	out = reflist_new();
Thomas White's avatar
Thomas White committed
176
	m = new_symopmask(holo);
Thomas White's avatar
Thomas White committed
177
178
179
180

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

182
		double total, sigma;
183
184
185
186
		signed int h, k, l;
		int n, j;
		int skip;

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

189
		/* There is a many-to-one correspondence between reflections
Thomas White's avatar
Thomas White committed
190
191
192
		 * in the merohedral and holohedral groups.  Do the calculation
		 * only once for each reflection in the holohedral group, which
		 * contains fewer reflections.
193
		 */
Thomas White's avatar
Thomas White committed
194
		get_asymm(holo, h, k, l, &h, &k, &l);
Thomas White's avatar
Thomas White committed
195
		if ( find_refl(out, h, k, l) != NULL ) continue;
196

197
198
		total = 0.0;
		sigma = 0.0;
199
		skip = 0;
Thomas White's avatar
Thomas White committed
200
201
202
		special_position(holo, m, h, k, l);
		n = num_equivs(holo, m);

203
204
205
206
207
		for ( j=0; j<n; j++ ) {

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

Thomas White's avatar
Thomas White committed
208
			get_equiv(holo, m, j, h, k, l, &he, &ke, &le);
209
210
211
212
213
214

			/* 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
215
216
			if ( !find_equiv_in_list(in, he, ke, le, mero,
			                         &hu, &ku, &lu) ) {
217
218
219
220
221
				/* 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
222
223
				      h, k, l, he, ke, le, symmetry_name(mero),
				      h, k, l);
224
225
226
227
				skip = 1;
				break;
			}

Thomas White's avatar
Thomas White committed
228
229
			total += get_intensity(refl);
			sigma += pow(get_esd_intensity(refl), 2.0);
230
231
232
233
234

		}

		if ( !skip ) {

Thomas White's avatar
Thomas White committed
235
236
237
			Reflection *new = add_refl(out, h, k, l);
			set_int(new, total);
			set_esd_intensity(new, sqrt(sigma));
238
239
240
241
242

		}

	}

Thomas White's avatar
Thomas White committed
243
	return out;
244
245
246
}


Thomas White's avatar
Thomas White committed
247
248
static RefList *expand_reflections(RefList *in, const SymOpList *target,
                                                const SymOpList *initial)
Thomas White's avatar
Thomas White committed
249
{
Thomas White's avatar
Thomas White committed
250
251
252
	Reflection *refl;
	RefListIterator *iter;
	RefList *out;
Thomas White's avatar
Thomas White committed
253
	SymOpMask *m;
Thomas White's avatar
Thomas White committed
254

Thomas White's avatar
Thomas White committed
255
	/* FIXME: Check properly */
Thomas White's avatar
Thomas White committed
256
	if ( num_equivs(target, NULL) > num_equivs(initial, NULL) ) {
Thomas White's avatar
Thomas White committed
257
258
		ERROR("%s is not a subgroup of %s!\n", symmetry_name(initial),
		                                       symmetry_name(target));
259
260
261
		return NULL;
	}

Thomas White's avatar
Thomas White committed
262
	out = reflist_new();
Thomas White's avatar
Thomas White committed
263
	m = new_symopmask(initial);
Thomas White's avatar
Thomas White committed
264
265
266
267

	for ( refl = first_refl(in, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
Thomas White's avatar
Thomas White committed
268
269
270
271
272

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

Thomas White's avatar
Thomas White committed
273
274
		get_indices(refl, &h, &k, &l);
		intensity = get_intensity(refl);
Thomas White's avatar
Thomas White committed
275

Thomas White's avatar
Thomas White committed
276
277
		special_position(initial, m, h, k, l);
		n = num_equivs(initial, m);
Thomas White's avatar
Thomas White committed
278

Thomas White's avatar
Thomas White committed
279
		/* For each equivalent in the higher symmetry group */
Thomas White's avatar
Thomas White committed
280
281
282
		for ( j=0; j<n; j++ ) {

			signed int he, ke, le;
Thomas White's avatar
Thomas White committed
283
			Reflection *new;
284
285
			int have_phase;
			double ph;
Thomas White's avatar
Thomas White committed
286
287

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

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

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

Thomas White's avatar
Thomas White committed
297
298
			/* FIXME: Make phase negative if the reflection is
			 * separated from the original via an inversion */
299
300
301
302
			get_phase(refl, &have_phase);
			if ( have_phase ) {
				set_ph(new, -ph);
			}
Thomas White's avatar
Thomas White committed
303

Thomas White's avatar
Thomas White committed
304
305
306
307
		}

	}

Thomas White's avatar
Thomas White committed
308
309
	free_symopmask(m);

Thomas White's avatar
Thomas White committed
310
	return out;
Thomas White's avatar
Thomas White committed
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
338
339
340
341
342
343
344
345
346
347
348
349
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;
}


350
351
352
int main(int argc, char *argv[])
{
	int c;
353
354
	int config_noise = 0;
	int config_poisson = 0;
355
	int config_multi = 0;
356
	int config_trimc = 0;
Thomas White's avatar
Thomas White committed
357
358
359
360
361
362
	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
363
364
	char *input_file = NULL;
	char *template = NULL;
Thomas White's avatar
Thomas White committed
365
	char *output = NULL;
366
	char *beamfile = NULL;
367
	struct beam_params *beam = NULL;
Thomas White's avatar
Thomas White committed
368
369
	RefList *input;
	UnitCell *cell = NULL;
370
371
372
373

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
Thomas White's avatar
Thomas White committed
374
		{"template",           1, NULL,               't'},
375
376
		{"poisson",            0, &config_poisson,     1},
		{"noise",              0, &config_noise,       1},
Thomas White's avatar
Thomas White committed
377
		{"output",             1, NULL,               'o'},
378
		{"symmetry",           1, NULL,               'y'},
Thomas White's avatar
Thomas White committed
379
380
		{"twin",               1, NULL,               'w'},
		{"expand",             1, NULL,               'e'},
381
		{"intensities",        1, NULL,               'i'},
382
		{"multiplicity",       0, &config_multi,       1},
383
		{"beam",               1, NULL,               'b'},
Thomas White's avatar
Thomas White committed
384
		{"pdb",                1, NULL,               'p'},
385
		{"trim-centrics",      0, &config_trimc,       1},
386
387
388
389
		{0, 0, NULL, 0}
	};

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

		switch (c) {
Thomas White's avatar
Thomas White committed
394
		case 'h' :
395
396
397
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
398
		case 't' :
Thomas White's avatar
Thomas White committed
399
400
401
			template = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
402
		case 'o' :
Thomas White's avatar
Thomas White committed
403
404
405
			output = strdup(optarg);
			break;

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

410
		case 'y' :
Thomas White's avatar
Thomas White committed
411
			mero_str = strdup(optarg);
412
413
414
			break;

		case 'w' :
Thomas White's avatar
Thomas White committed
415
			holo_str = strdup(optarg);
416
417
			break;

Thomas White's avatar
Thomas White committed
418
		case 'e' :
Thomas White's avatar
Thomas White committed
419
			expand_str = strdup(optarg);
Thomas White's avatar
Thomas White committed
420
421
			break;

422
423
424
425
		case 'b' :
			beamfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
426
427
428
429
430
431
432
433
		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
434
		case 0 :
435
436
			break;

Thomas White's avatar
Thomas White committed
437
		default :
438
439
440
441
442
			return 1;
		}

	}

Thomas White's avatar
Thomas White committed
443
	if ( (holo_str != NULL) && (expand_str != NULL) ) {
Thomas White's avatar
Thomas White committed
444
445
		ERROR("You cannot 'twin' and 'expand' at the same time.\n");
		ERROR("Decide which one you want to do first.\n");
446
447
448
449
450
451
452
453
454
455
		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
456
457
	}

Thomas White's avatar
Thomas White committed
458
459
460
	if ( cell == NULL ) {
		ERROR("You need to give a PDB file with the unit cell.\n");
		return 1;
461
	}
462

Thomas White's avatar
Thomas White committed
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
	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
482
	input = read_reflections(input_file);
Thomas White's avatar
Thomas White committed
483
484
485
486
	if ( input == NULL ) {
		ERROR("Problem reading input file %s\n", input_file);
		return 1;
	}
Thomas White's avatar
Thomas White committed
487
	free(input_file);
488
489
490
491

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

	if ( !config_trimc && check_list_symmetry(input, mero) ) {
492
		ERROR("The input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
493
		      " have symmetry %s\n", symmetry_name(mero));
494
		return 1;
495
	}
496

497
498
	if ( config_poisson ) {
		if ( beam != NULL ) {
Thomas White's avatar
Thomas White committed
499
			poisson_reflections(input, beam->adu_per_photon);
500
501
502
503
504
505
506
		} 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
507
	if ( config_noise ) noise_reflections(input);
Thomas White's avatar
Thomas White committed
508

509
	if ( holo != NULL ) {
510

Thomas White's avatar
Thomas White committed
511
		RefList *new;
Thomas White's avatar
Thomas White committed
512
513
		STATUS("Twinning from %s into %s\n", symmetry_name(mero),
		                                     symmetry_name(holo));
Thomas White's avatar
Thomas White committed
514
515
516
517
518
519
520
521
522
		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;
523

Thomas White's avatar
Thomas White committed
524
525
	}

Thomas White's avatar
Thomas White committed
526
	if ( expand != NULL ) {
527

Thomas White's avatar
Thomas White committed
528
		RefList *new;
Thomas White's avatar
Thomas White committed
529
530
		STATUS("Expanding from %s into %s\n", symmetry_name(mero),
		                                      symmetry_name(expand));
Thomas White's avatar
Thomas White committed
531
532
533
534
535
		new = expand_reflections(input, expand, mero);

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

Thomas White's avatar
Thomas White committed
537
538
	}

539
540
541
542
543
544
545
546
547
548
549
550
551
	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));

	}

552
553
	if ( config_multi ) {

Thomas White's avatar
Thomas White committed
554
555
		Reflection *refl;
		RefListIterator *iter;
Thomas White's avatar
Thomas White committed
556
557
558
		SymOpMask *m;

		m = new_symopmask(mero);
559

Thomas White's avatar
Thomas White committed
560
561
562
		for ( refl = first_refl(input, &iter);
		      refl != NULL;
		      refl = next_refl(refl, iter) ) {
563
564

			double inty;
Thomas White's avatar
Thomas White committed
565
566
567
568
			signed int h, k, l;

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

Thomas White's avatar
Thomas White committed
570
571
			special_position(mero, m, h, k, l);
			inty *= (double)num_equivs(mero, m);
Thomas White's avatar
Thomas White committed
572
			set_int(refl, inty);
573
574

		}
Thomas White's avatar
Thomas White committed
575
576

		free_symopmask(m);
577
578
	}

Thomas White's avatar
Thomas White committed
579
	if ( template ) {
580

Thomas White's avatar
Thomas White committed
581
582
583
584
		RefList *t = read_reflections(template);
		RefList *new = template_reflections(input, t);
		reflist_free(input);
		input = new;
585

Thomas White's avatar
Thomas White committed
586
587
	}

Thomas White's avatar
Thomas White committed
588
	write_reflist(output, input, cell);
Thomas White's avatar
Thomas White committed
589

Thomas White's avatar
Thomas White committed
590
	reflist_free(input);
591
592
593

	return 0;
}