compare_hkl.c 11.5 KB
Newer Older
1
2
3
4
5
/*
 * compare_hkl.c
 *
 * Compare reflection lists
 *
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>
24
25

#include "utils.h"
Thomas White's avatar
Thomas White committed
26
#include "statistics.h"
27
#include "symmetry.h"
Thomas White's avatar
Thomas White committed
28
#include "reflist-utils.h"
29
30


Thomas White's avatar
Thomas White committed
31
/* Number of bins for plot of resolution shells */
32
#define NBINS (10)
Thomas White's avatar
Thomas White committed
33
34


35
36
static void show_help(const char *s)
{
37
	printf("Syntax: %s [options] <file1.hkl> <file2.hkl>\n\n", s);
38
39
40
41
	printf(
"Compare intensity lists.\n"
"\n"
"  -h, --help                 Display this help message.\n"
Thomas White's avatar
Thomas White committed
42
"  -o, --ratio=<filename>     Specify output filename for ratios.\n"
43
"  -y, --symmetry=<sym>       The symmetry of both the input files.\n"
44
"  -p, --pdb=<filename>       PDB file to use (default: molecule.pdb).\n"
Thomas White's avatar
Thomas White committed
45
"      --shells               Plot the figures of merit by resolution.\n"
46
47
"      --rmin=<res>           Fix lower resolution limit for --shells (m^-1).\n"
"      --rmax=<res>           Fix upper resolution limit for --shells (m^-1).\n"
48
49
50
51
"\n");
}


Thomas White's avatar
Thomas White committed
52
53
54
static void plot_shells(RefList *list1, RefList *list2, double scale,
                        UnitCell *cell, const char *sym,
                        double rmin_fix, double rmax_fix)
Thomas White's avatar
Thomas White committed
55
56
{
	double num[NBINS];
Thomas White's avatar
Thomas White committed
57
	int cts[NBINS];
58
59
60
61
62
63
64
	int possible[NBINS];
	unsigned int *counted;
	unsigned int measurements[NBINS];
	unsigned int measured[NBINS];
	double total_vol, vol_per_shell;
	double rmins[NBINS];
	double rmaxs[NBINS];
65
	double snr[NBINS];
Thomas White's avatar
Thomas White committed
66
	double den;
67
68
	double rmin, rmax;
	signed int h, k, l;
Thomas White's avatar
Thomas White committed
69
	int i;
Thomas White's avatar
Thomas White committed
70
	int ctot;
Thomas White's avatar
Thomas White committed
71
	FILE *fh;
72
	int nout = 0;
Thomas White's avatar
Thomas White committed
73
74
75
76
77
78
	Reflection *refl1;
	RefListIterator *iter;
	double asx, asy, asz;
	double bsx, bsy, bsz;
	double csx, csy, csz;
	int hmax, kmax, lmax;
Thomas White's avatar
Thomas White committed
79
80

	if ( cell == NULL ) {
Thomas White's avatar
Thomas White committed
81
		ERROR("Need the unit cell to plot resolution shells.\n");
Thomas White's avatar
Thomas White committed
82
83
84
		return;
	}

Thomas White's avatar
Thomas White committed
85
	fh = fopen("shells.dat", "w");
Thomas White's avatar
Thomas White committed
86
	if ( fh == NULL ) {
Thomas White's avatar
Thomas White committed
87
		ERROR("Couldn't open 'shells.dat'\n");
Thomas White's avatar
Thomas White committed
88
89
90
91
92
		return;
	}

	for ( i=0; i<NBINS; i++ ) {
		num[i] = 0.0;
Thomas White's avatar
Thomas White committed
93
		cts[i] = 0;
94
95
96
		possible[i] = 0;
		measured[i] = 0;
		measurements[i] = 0;
97
		snr[i] = 0;
Thomas White's avatar
Thomas White committed
98
99
	}

Thomas White's avatar
Thomas White committed
100
101
102
103
104
105
	/* Iterate over all common reflections and calculate min and max
	 * resolution */
	rmin = +INFINITY;  rmax = 0.0;
	for ( refl1 = first_refl(list1, &iter);
	      refl1 != NULL;
	      refl1 = next_refl(refl1, iter) ) {
Thomas White's avatar
Thomas White committed
106
107

		signed int h, k, l;
108
		double d;
Thomas White's avatar
Thomas White committed
109
		Reflection *refl2;
Thomas White's avatar
Thomas White committed
110

Thomas White's avatar
Thomas White committed
111
112
113
		get_indices(refl1, &h, &k, &l);
		refl2 = find_refl(list2, h, k, l);
		if ( refl2 == NULL ) continue;  /* No common reflection */
Thomas White's avatar
Thomas White committed
114

115
		d = resolution(cell, h, k, l) * 2.0;
116
117
		if ( d > rmax ) rmax = d;
		if ( d < rmin ) rmin = d;
Thomas White's avatar
Thomas White committed
118
119

	}
120

121
	STATUS("1/d goes from %f to %f nm^-1\n", rmin/1e9, rmax/1e9);
122

123
124
	/* Widen the range just a little bit */
	rmin -= 0.001e9;
125
126
	rmax += 0.001e9;

127
128
129
	/* Fixed resolution shells if needed */
	if ( rmin_fix > 0.0 ) rmin = rmin_fix;
	if ( rmax_fix > 0.0 ) rmax = rmax_fix;
130

Thomas White's avatar
Thomas White committed
131
	/* Calculate the resolution bins */
132
133
134
135
136
137
138
139
140
141
	total_vol = pow(rmax, 3.0) - pow(rmin, 3.0);
	vol_per_shell = total_vol / NBINS;
	rmins[0] = rmin;
	for ( i=1; i<NBINS; i++ ) {

		double r;

		r = vol_per_shell + pow(rmins[i-1], 3.0);
		r = pow(r, 1.0/3.0);

142
		/* Shells of constant volume */
143
144
145
		rmaxs[i-1] = r;
		rmins[i] = r;

146
147
148
149
		/* Shells of constant thickness */
		//rmins[i] = rmins[i-1] + (rmax-rmin)/NBINS;
		//rmaxs[i-1] = rmins[i-1] + (rmax-rmin)/NBINS;

150
151
152
153
154
155
156
157
158
159
		STATUS("Shell %i: %f to %f\n", i-1,
		       rmins[i-1]/1e9, rmaxs[i-1]/1e9);

	}
	rmaxs[NBINS-1] = rmax;
	STATUS("Shell %i: %f to %f\n", NBINS-1,
	       rmins[NBINS-1]/1e9, rmaxs[NBINS-1]/1e9);

	/* Count the number of reflections possible in each shell */
	counted = new_list_count();
Thomas White's avatar
Thomas White committed
160
161
162
163
164
165
166
167
168
169
	cell_get_reciprocal(cell, &asx, &asy, &asz,
	                          &bsx, &bsy, &bsz,
	                          &csx, &csy, &csz);
	hmax = rmax / modulus(asx, asy, asz);
	kmax = rmax / modulus(bsx, bsy, bsz);
	lmax = rmax / modulus(csx, csy, csz);

	for ( h=-hmax; h<hmax; h++ ) {
	for ( k=-kmax; k<kmax; k++ ) {
	for ( l=-lmax; l<lmax; l++ ) {
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196

		double d;
		signed int hs, ks, ls;
		int bin;

		get_asymm(h, k, l, &hs, &ks, &ls, sym);
		if ( lookup_count(counted, hs, ks, ls) ) continue;
		set_count(counted, hs, ks, ls, 1);

		d = resolution(cell, h, k, l) * 2.0;

		bin = -1;
		for ( i=0; i<NBINS; i++ ) {
			if ( (d>rmins[i]) && (d<=rmaxs[i]) ) {
				bin = i;
				break;
			}
		}
		if ( bin == -1 ) continue;

		possible[bin]++;

	}
	}
	}
	free(counted);

Thomas White's avatar
Thomas White committed
197
	den = 0.0;
Thomas White's avatar
Thomas White committed
198
	ctot = 0;
199
	nout = 0;
Thomas White's avatar
Thomas White committed
200

Thomas White's avatar
Thomas White committed
201
202
203
204
	for ( refl1 = first_refl(list1, &iter);
	      refl1 != NULL;
	      refl1 = next_refl(refl1, iter) ) {

Thomas White's avatar
Thomas White committed
205
		signed int h, k, l;
206
		double d;
Thomas White's avatar
Thomas White committed
207
		int bin;
208
		double i1, i2;
209
		int j;
Thomas White's avatar
Thomas White committed
210
		Reflection *refl2;
Thomas White's avatar
Thomas White committed
211

Thomas White's avatar
Thomas White committed
212
213
214
		get_indices(refl1, &h, &k, &l);
		refl2 = find_refl(list2, h, k, l);
		if ( refl2 == NULL ) continue;  /* No common reflection */
Thomas White's avatar
Thomas White committed
215

216
		d = resolution(cell, h, k, l) * 2.0;
Thomas White's avatar
Thomas White committed
217

218
219
220
221
222
223
224
		bin = -1;
		for ( j=0; j<NBINS; j++ ) {
			if ( (d>rmins[j]) && (d<=rmaxs[j]) ) {
				bin = j;
				break;
			}
		}
225
226

		/* Outside resolution range? */
227
228
229
230
		if ( bin == -1 ) {
			nout++;
			continue;
		}
Thomas White's avatar
Thomas White committed
231

Thomas White's avatar
Thomas White committed
232
233
		i1 = get_intensity(refl1);
		i2 = get_intensity(refl2);
Thomas White's avatar
Thomas White committed
234
		i2 *= scale;
Thomas White's avatar
Thomas White committed
235

Thomas White's avatar
Thomas White committed
236
		num[bin] += fabs(i1 - i2);
237
		den += i1;
Thomas White's avatar
Thomas White committed
238
		ctot++;
239
		cts[bin]++;
Thomas White's avatar
Thomas White committed
240
241
242

	}

243
244
245
246
247
	if ( nout ) {
		STATUS("Warning; %i reflections outside resolution range.\n",
		       nout);
	}

Thomas White's avatar
Thomas White committed
248
249
	for ( i=0; i<NBINS; i++ ) {

250
		double r, cen;
251
		cen = rmins[i] + (rmaxs[i] - rmins[i])/2.0;
252
		r = (num[i]/den)*((double)ctot/cts[i]);
253
		fprintf(fh, "%f %f\n", cen*1.0e-9, r*100.0);
Thomas White's avatar
Thomas White committed
254
255
256
257
258
259
260

	}

	fclose(fh);
}


261
262
263
int main(int argc, char *argv[])
{
	int c;
264
	UnitCell *cell;
Thomas White's avatar
Thomas White committed
265
	char *ratiofile = NULL;
266
267
	char *afile = NULL;
	char *bfile = NULL;
Thomas White's avatar
Thomas White committed
268
	char *sym = NULL;
269
	double scale, scale_r2, scale_rdig, R1, R2, R1i, Rdiff, pearson;
270
	double scale_rintint, scale_r1i, scale_r1, scale_r1fi;
Thomas White's avatar
Thomas White committed
271
272
273
274
275
276
	int ncom;
	RefList *list1;
	RefList *list2;
	RefList *list2_transformed;
	RefList *ratio;
	RefList *deleteme;
Thomas White's avatar
Thomas White committed
277
	int config_shells = 0;
278
	char *pdb = NULL;
279
280
	int rej1 = 0;
	int rej2 = 0;
281
282
	float rmin_fix = -1.0;
	float rmax_fix = -1.0;
Thomas White's avatar
Thomas White committed
283
284
	Reflection *refl1;
	RefListIterator *iter;
285
286
287
288

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
Thomas White's avatar
Thomas White committed
289
		{"ratio" ,             1, NULL,               'o'},
290
		{"symmetry",           1, NULL,               'y'},
Thomas White's avatar
Thomas White committed
291
		{"shells",             0, &config_shells,     1},
292
		{"pdb",                1, NULL,               'p'},
293
294
		{"rmin",               1, NULL,               2},
		{"rmax",               1, NULL,               3},
295
296
297
298
		{0, 0, NULL, 0}
	};

	/* Short options */
299
	while ((c = getopt_long(argc, argv, "ho:y:p:", longopts, NULL)) != -1) {
300
301

		switch (c) {
Thomas White's avatar
Thomas White committed
302
		case 'h' :
303
304
305
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
306
		case 'o' :
Thomas White's avatar
Thomas White committed
307
			ratiofile = strdup(optarg);
308
309
			break;

310
311
312
313
		case 'y' :
			sym = strdup(optarg);
			break;

314
315
316
317
		case 'p' :
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
318
		case 0 :
319
320
			break;

321
322
323
324
325
326
327
328
329
330
331
332
333
334
		case 2 :
			if ( sscanf(optarg, "%e", &rmin_fix) != 1 ) {
				ERROR("Invalid value for --rmin\n");
				return 1;
			}
			break;

		case 3 :
			if ( sscanf(optarg, "%e", &rmax_fix) != 1 ) {
				ERROR("Invalid value for --rmax\n");
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
335
		default :
336
337
338
339
340
			return 1;
		}

	}

341
342
343
344
345
	if ( argc != (optind+2) ) {
		ERROR("Please provide exactly two HKL files to compare.\n");
		return 1;
	}

346
347
348
349
	if ( sym == NULL ) {
		sym = strdup("1");
	}

350
351
352
	afile = strdup(argv[optind++]);
	bfile = strdup(argv[optind]);

353
354
355
356
357
358
359
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

	cell = load_cell_from_pdb(pdb);
	free(pdb);

Thomas White's avatar
Thomas White committed
360
361
362
	list1 = read_reflections(afile);
	if ( list1 == NULL ) {
		ERROR("Couldn't read file '%s'\n", afile);
363
364
		return 1;
	}
Thomas White's avatar
Thomas White committed
365
366
367
368

	list2 = read_reflections(bfile);
	if ( list2 == NULL ) {
		ERROR("Couldn't read file '%s'\n", bfile);
369
370
371
		return 1;
	}

Thomas White's avatar
Thomas White committed
372
	/* Check that the intensities have the correct symmetry */
Thomas White's avatar
Thomas White committed
373
	if ( check_list_symmetry(list1, sym) ) {
Thomas White's avatar
Thomas White committed
374
375
376
377
		ERROR("The first input reflection list does not appear to"
		      " have symmetry %s\n", sym);
		return 1;
	}
Thomas White's avatar
Thomas White committed
378
	if ( check_list_symmetry(list2, sym) ) {
Thomas White's avatar
Thomas White committed
379
380
381
382
383
		ERROR("The second input reflection list does not appear to"
		      " have symmetry %s\n", sym);
		return 1;
	}

384
	/* Find common reflections (taking symmetry into account) */
Thomas White's avatar
Thomas White committed
385
386
387
388
389
390
	list2_transformed = reflist_new();
	ratio = reflist_new();
	deleteme = reflist_new();
	for ( refl1 = first_refl(list1, &iter);
	      refl1 != NULL;
	      refl1 = next_refl(refl1, iter) ) {
391

Thomas White's avatar
Thomas White committed
392
		signed int h, k, l;
393
394
		signed int he, ke, le;
		double val1, val2;
395
396
		double sig1, sig2;
		int ig = 0;
397
		double d;
Thomas White's avatar
Thomas White committed
398
399
		Reflection *refl2;
		Reflection *tr;
400

Thomas White's avatar
Thomas White committed
401
		get_indices(refl1, &h, &k, &l);
402

Thomas White's avatar
Thomas White committed
403
404
405
		if ( !find_equiv_in_list(list2, h, k, l, sym, &he, &ke, &le) ) {
			/* No common reflection */
			add_refl(deleteme, h, k, l);
406
407
			continue;
		}
408

Thomas White's avatar
Thomas White committed
409
410
411
412
413
414
415
		refl2 = find_refl(list2, he, ke, le);
		assert(refl2 != NULL);

		val1 = get_intensity(refl1);
		val2 = get_intensity(refl2);
		sig1 = get_esd_intensity(refl1);
		sig2 = get_esd_intensity(refl2);
416
417
418
419
420
421
422
423
424

		if ( val1 < 3.0 * sig1 ) {
			rej1++;
			ig = 1;
		}
		if ( val2 < 3.0 * sig2 ) {
			rej2++;
			ig = 1;
		}
425
426
427
428
429

		d = 0.5/resolution(cell, h, k, l);
		if ( d > 55.0e-10 ) ig = 1;
		//if ( d < 15.0e-10 ) ig = 1;

430
		//if ( ig ) continue;
431

Thomas White's avatar
Thomas White committed
432
433
434
435
436
437
438
439
		/* Add the old data from 'refl2' to a new list with the same
		 * indices as its equivalent in 'list1' */
		tr = add_refl(list2_transformed, h, k, l);
		copy_data(tr, refl2);

		/* Add divided version to 'output' list */
		tr = add_refl(ratio, h, k, l);
		set_int(tr, val1/val2);
440
441

	}
Thomas White's avatar
Thomas White committed
442

Thomas White's avatar
Thomas White committed
443
444
	STATUS("%i reflections in '%s' had I < 3.0*sigma(I)\n", rej1, afile);
	STATUS("%i reflections in '%s' had I < 3.0*sigma(I)\n", rej2, bfile);
445

Thomas White's avatar
Thomas White committed
446
	ncom = num_reflections(list2_transformed);
Thomas White's avatar
Thomas White committed
447
	STATUS("%i,%i reflections: %i in common\n",
Thomas White's avatar
Thomas White committed
448
449
450
451
452
453
	       num_reflections(list1), num_reflections(list2), ncom);

	/* Trim reflections from 'list1' which had no equivalents in 'list2' */
	for ( refl1 = first_refl(deleteme, &iter);
	      refl1 != NULL;
	      refl1 = next_refl(refl1, iter) ) {
Thomas White's avatar
Thomas White committed
454

Thomas White's avatar
Thomas White committed
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
		signed int h, k, l;
		Reflection *del;

		get_indices(refl1, &h, &k, &l);
		del = find_refl(list1, h, k, l);
		assert(del != NULL);

		delete_refl(del);

	}
	reflist_free(deleteme);
	reflist_free(list2);

	R1 = stat_r1_ignore(list1, list2_transformed, &scale_r1fi);
	STATUS("R1(F) = %5.4f %% (scale=%5.2e)"
	       " (ignoring negative intensities)\n",
471
	       R1*100.0, scale_r1fi);
Thomas White's avatar
Thomas White committed
472

Thomas White's avatar
Thomas White committed
473
474
475
	R1 = stat_r1_zero(list1, list2_transformed, &scale_r1);
	STATUS("R1(F) = %5.4f %% (scale=%5.2e)"
	       " (zeroing negative intensities)\n",
476
	       R1*100.0, scale_r1);
Thomas White's avatar
Thomas White committed
477

Thomas White's avatar
Thomas White committed
478
	R2 = stat_r2(list1, list2_transformed, &scale_r2);
479
	STATUS("R2(I) = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale_r2);
Thomas White's avatar
Thomas White committed
480

Thomas White's avatar
Thomas White committed
481
	R1i = stat_r1_i(list1, list2_transformed, &scale_r1i);
482
	STATUS("R1(I) = %5.4f %% (scale=%5.2e)\n", R1i*100.0, scale_r1i);
Thomas White's avatar
Thomas White committed
483

Thomas White's avatar
Thomas White committed
484
485
486
	Rdiff = stat_rdiff_ignore(list1, list2_transformed, &scale_rdig);
	STATUS("Rint(F) = %5.4f %% (scale=%5.2e)"
	       " (ignoring negative intensities)\n",
487
	       Rdiff*100.0, scale_rdig);
Thomas White's avatar
Thomas White committed
488

Thomas White's avatar
Thomas White committed
489
490
491
	Rdiff = stat_rdiff_zero(list1, list2_transformed, &scale);
	STATUS("Rint(F) = %5.4f %% (scale=%5.2e)"
	       " (zeroing negative intensities)\n",
Thomas White's avatar
Thomas White committed
492
493
	       Rdiff*100.0, scale);

Thomas White's avatar
Thomas White committed
494
	Rdiff = stat_rdiff_intensity(list1, list2_transformed, &scale_rintint);
Thomas White's avatar
Thomas White committed
495
	STATUS("Rint(I) = %5.4f %% (scale=%5.2e)\n",
496
	       Rdiff*100.0, scale_rintint);
Thomas White's avatar
Thomas White committed
497

Thomas White's avatar
Thomas White committed
498
	pearson = stat_pearson_i(list1, list2_transformed);
Thomas White's avatar
Thomas White committed
499
500
	STATUS("Pearson r(I) = %5.4f\n", pearson);

Thomas White's avatar
Thomas White committed
501
	pearson = stat_pearson_f_ignore(list1, list2_transformed);
Thomas White's avatar
Thomas White committed
502
503
504
	STATUS("Pearson r(F) = %5.4f (ignoring negative intensities)\n",
	       pearson);

Thomas White's avatar
Thomas White committed
505
	pearson = stat_pearson_f_zero(list1, list2_transformed);
Thomas White's avatar
Thomas White committed
506
507
	STATUS("Pearson r(F) = %5.4f (zeroing negative intensities)\n",
	       pearson);
Thomas White's avatar
Thomas White committed
508

Thomas White's avatar
Thomas White committed
509
	if ( config_shells ) {
Thomas White's avatar
Thomas White committed
510
		plot_shells(list1, list2_transformed, scale_r1fi,
511
		            cell, sym, rmin_fix, rmax_fix);
Thomas White's avatar
Thomas White committed
512
513
	}

Thomas White's avatar
Thomas White committed
514
515
	if ( ratiofile != NULL ) {
		write_reflist(ratiofile, ratio, cell);
516
	}
517
518
519

	return 0;
}