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"
Thomas White's avatar
Thomas White committed
61
62
"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"
63
64
"  -w, --twin=<sym>           Generate twinned data according to the given\n"
"                              point group.\n"
Thomas White's avatar
Thomas White committed
65
"  -e, --expand=<sym>         Expand reflections to this point group.\n"
66
"\n"
67
68
69
70
71
"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"
72
73
74
75
"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"
76
77
"      --multiplicity         Multiply intensities by the number of\n"
"                              equivalent reflections.\n"
78
79
80
"\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
81
82
"  -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
83
);
Thomas White's avatar
Thomas White committed
84
85
86
87
}


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

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

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

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

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

	}
105
106
107
108
}


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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

		}

		if ( !skip ) {

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

		}

	}

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

	}

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

Thomas White's avatar
Thomas White committed
311
	return out;
Thomas White's avatar
Thomas White committed
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
350
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;
}


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

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

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

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

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

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

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

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

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

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

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

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

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

	}

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

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

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

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

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

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

510
	if ( holo != NULL ) {
511

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

Thomas White's avatar
Thomas White committed
525
526
	}

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

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

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

Thomas White's avatar
Thomas White committed
538
539
	}

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

	}

553
554
	if ( config_multi ) {

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

		m = new_symopmask(mero);
560

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

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

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

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

		}
Thomas White's avatar
Thomas White committed
576
577

		free_symopmask(m);
578
579
	}

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

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

Thomas White's avatar
Thomas White committed
587
588
	}

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

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

	return 0;
}