compare_hkl.c 10.6 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
/*
 * compare_hkl.c
 *
 * Compare reflection lists
 *
 * (c) 2006-2010 Thomas White <taw@physics.org>
 *
 * 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 "reflections.h"
Thomas White's avatar
Thomas White committed
26
#include "statistics.h"
27
#include "symmetry.h"
28
29


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


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


Thomas White's avatar
Thomas White committed
51
static void plot_shells(const double *ref1, const double *ref2,
52
                        ReflItemList *items, double scale, UnitCell *cell,
53
                        const char *sym, double rmin_fix, double rmax_fix)
Thomas White's avatar
Thomas White committed
54
55
{
	double num[NBINS];
Thomas White's avatar
Thomas White committed
56
	int cts[NBINS];
57
58
59
60
61
62
63
	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];
64
	double snr[NBINS];
Thomas White's avatar
Thomas White committed
65
	double den;
66
67
	double rmin, rmax;
	signed int h, k, l;
Thomas White's avatar
Thomas White committed
68
	int i;
Thomas White's avatar
Thomas White committed
69
	int ctot;
Thomas White's avatar
Thomas White committed
70
	FILE *fh;
71
	int nout = 0;
Thomas White's avatar
Thomas White committed
72
73

	if ( cell == NULL ) {
Thomas White's avatar
Thomas White committed
74
		ERROR("Need the unit cell to plot resolution shells.\n");
Thomas White's avatar
Thomas White committed
75
76
77
		return;
	}

Thomas White's avatar
Thomas White committed
78
	fh = fopen("shells.dat", "w");
Thomas White's avatar
Thomas White committed
79
	if ( fh == NULL ) {
Thomas White's avatar
Thomas White committed
80
		ERROR("Couldn't open 'shells.dat'\n");
Thomas White's avatar
Thomas White committed
81
82
83
84
85
		return;
	}

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

	rmin = +INFINITY;
	rmax = 0.0;
	for ( i=0; i<num_items(items); i++ ) {

		struct refl_item *it;
		signed int h, k, l;
99
		double d;
Thomas White's avatar
Thomas White committed
100
101
102
103

		it = get_item(items, i);
		h = it->h;  k = it->k;  l = it->l;

104
		d = resolution(cell, h, k, l) * 2.0;
105
106
		if ( d > rmax ) rmax = d;
		if ( d < rmin ) rmin = d;
Thomas White's avatar
Thomas White committed
107
108

	}
109

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

112
113
	/* Widen the range just a little bit */
	rmin -= 0.001e9;
114
115
	rmax += 0.001e9;

116
117
118
	/* Fixed resolution shells if needed */
	if ( rmin_fix > 0.0 ) rmin = rmin_fix;
	if ( rmax_fix > 0.0 ) rmax = rmax_fix;
119

120
121
122
123
124
125
126
127
128
129
	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);

130
		/* Shells of constant volume */
131
132
133
		rmaxs[i-1] = r;
		rmins[i] = r;

134
135
136
137
		/* Shells of constant thickness */
		//rmins[i] = rmins[i-1] + (rmax-rmin)/NBINS;
		//rmaxs[i-1] = rmins[i-1] + (rmax-rmin)/NBINS;

138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
		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();
	for ( h=-50; h<=+50; h++ ) {
	for ( k=-50; k<=+50; k++ ) {
	for ( l=-50; l<=+50; l++ ) {

		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
178
	den = 0.0;
Thomas White's avatar
Thomas White committed
179
	ctot = 0;
180
	nout = 0;
Thomas White's avatar
Thomas White committed
181
182
183
184
	for ( i=0; i<num_items(items); i++ ) {

		struct refl_item *it;
		signed int h, k, l;
185
		double d;
Thomas White's avatar
Thomas White committed
186
		int bin;
187
		double i1, i2;
188
		int j;
Thomas White's avatar
Thomas White committed
189
190
191
192

		it = get_item(items, i);
		h = it->h;  k = it->k;  l = it->l;

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

195
196
197
198
199
200
201
		bin = -1;
		for ( j=0; j<NBINS; j++ ) {
			if ( (d>rmins[j]) && (d<=rmaxs[j]) ) {
				bin = j;
				break;
			}
		}
202
203

		/* Outside resolution range? */
204
205
206
207
		if ( bin == -1 ) {
			nout++;
			continue;
		}
Thomas White's avatar
Thomas White committed
208
209

		i1 = lookup_intensity(ref1, h, k, l);
210
		i2 = lookup_intensity(ref2, h, k, l);
Thomas White's avatar
Thomas White committed
211
		i2 *= scale;
Thomas White's avatar
Thomas White committed
212

Thomas White's avatar
Thomas White committed
213
		num[bin] += fabs(i1 - i2);
214
		den += i1;
Thomas White's avatar
Thomas White committed
215
		ctot++;
216
		cts[bin]++;
Thomas White's avatar
Thomas White committed
217
218
219

	}

220
221
222
223
224
	if ( nout ) {
		STATUS("Warning; %i reflections outside resolution range.\n",
		       nout);
	}

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

227
		double r, cen;
228
		cen = rmins[i] + (rmaxs[i] - rmins[i])/2.0;
229
		r = (num[i]/den)*((double)ctot/cts[i]);
230
		fprintf(fh, "%f %f\n", cen*1.0e-9, r*100.0);
Thomas White's avatar
Thomas White committed
231
232
233
234
235
236
237

	}

	fclose(fh);
}


238
239
240
241
242
int main(int argc, char *argv[])
{
	int c;
	double *ref1;
	double *ref2;
243
	double *ref2_transformed;
244
	double *out;
245
	UnitCell *cell;
246
247
248
	char *outfile = NULL;
	char *afile = NULL;
	char *bfile = NULL;
Thomas White's avatar
Thomas White committed
249
	char *sym = NULL;
250
	double scale, scale_r2, scale_rdig, R1, R2, R1i, Rdiff, pearson;
251
	double scale_rintint, scale_r1i, scale_r1, scale_r1fi;
Thomas White's avatar
Thomas White committed
252
253
	int i, ncom;
	ReflItemList *i1, *i2, *icommon;
Thomas White's avatar
Thomas White committed
254
	int config_shells = 0;
255
	char *pdb = NULL;
256
257
258
259
	double *esd1;
	double *esd2;
	int rej1 = 0;
	int rej2 = 0;
260
	unsigned int *cts1;
261
262
	float rmin_fix = -1.0;
	float rmax_fix = -1.0;
263
264
265
266
267

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"output",             1, NULL,               'o'},
268
		{"symmetry",           1, NULL,               'y'},
Thomas White's avatar
Thomas White committed
269
		{"shells",             0, &config_shells,     1},
270
		{"pdb",                1, NULL,               'p'},
271
272
		{"rmin",               1, NULL,               2},
		{"rmax",               1, NULL,               3},
273
274
275
276
		{0, 0, NULL, 0}
	};

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

		switch (c) {
Thomas White's avatar
Thomas White committed
280
		case 'h' :
281
282
283
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
284
		case 'o' :
285
286
287
			outfile = strdup(optarg);
			break;

288
289
290
291
		case 'y' :
			sym = strdup(optarg);
			break;

292
293
294
295
		case 'p' :
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
296
		case 0 :
297
298
			break;

299
300
301
302
303
304
305
306
307
308
309
310
311
312
		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
313
		default :
314
315
316
317
318
			return 1;
		}

	}

319
320
321
322
323
	if ( argc != (optind+2) ) {
		ERROR("Please provide exactly two HKL files to compare.\n");
		return 1;
	}

324
325
326
327
	if ( sym == NULL ) {
		sym = strdup("1");
	}

328
329
330
	afile = strdup(argv[optind++]);
	bfile = strdup(argv[optind]);

331
332
333
334
335
336
337
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

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

Thomas White's avatar
Thomas White committed
338
	ref1 = new_list_intensity();
339
	esd1 = new_list_sigma();
340
341
	cts1 = new_list_count();
	i1 = read_reflections(afile, ref1, NULL, cts1, esd1);
Thomas White's avatar
Thomas White committed
342
	if ( i1 == NULL ) {
343
344
345
		ERROR("Couldn't open file '%s'\n", afile);
		return 1;
	}
Thomas White's avatar
Thomas White committed
346
	ref2 = new_list_intensity();
347
348
	esd2 = new_list_sigma();
	i2 = read_reflections(bfile, ref2, NULL, NULL, esd2);
Thomas White's avatar
Thomas White committed
349
	if ( i2 == NULL ) {
350
351
352
353
		ERROR("Couldn't open file '%s'\n", bfile);
		return 1;
	}

Thomas White's avatar
Thomas White committed
354
355
356
357
358
359
360
361
362
363
364
365
	/* Check that the intensities have the correct symmetry */
	if ( check_symmetry(i1, sym) ) {
		ERROR("The first input reflection list does not appear to"
		      " have symmetry %s\n", sym);
		return 1;
	}
	if ( check_symmetry(i2, sym) ) {
		ERROR("The second input reflection list does not appear to"
		      " have symmetry %s\n", sym);
		return 1;
	}

Thomas White's avatar
Thomas White committed
366
367
	/* List for output scale factor map */
	out = new_list_intensity();
368

369
370
371
372
373
	/* Find common reflections (taking symmetry into account) */
	icommon = new_items();
	ref2_transformed = new_list_intensity();
	for ( i=0; i<num_items(i1); i++ ) {

Thomas White's avatar
Thomas White committed
374
375
		struct refl_item *it;
		signed int h, k, l;
376
377
		signed int he, ke, le;
		double val1, val2;
378
379
		double sig1, sig2;
		int ig = 0;
380
		double d;
381

382
		it = get_item(i1, i);
Thomas White's avatar
Thomas White committed
383
		h = it->h;  k = it->k;  l = it->l;
384

385
		if ( !find_unique_equiv(i2, h, k, l, sym, &he, &ke, &le) ) {
Thomas White's avatar
Thomas White committed
386
387
			//STATUS("%i %i %i not matched (%f nm).\n", h, k, l,
			//       1.0/(2.0*resolution(cell, h, k, l)/1e9));
388
389
			continue;
		}
390

391
392
		val1 = lookup_intensity(ref1, h, k, l);
		val2 = lookup_intensity(ref2, he, ke, le);
393
		sig1 = lookup_sigma(esd1, h, k, l);
394
		sig2 = lookup_sigma(esd2, he, ke, le);
395
396
397
398
399
400
401
402
403

		if ( val1 < 3.0 * sig1 ) {
			rej1++;
			ig = 1;
		}
		if ( val2 < 3.0 * sig2 ) {
			rej2++;
			ig = 1;
		}
404
405
406
407
408

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

409
		//if ( ig ) continue;
410

411
412
413
		set_intensity(ref2_transformed, h, k, l, val2);
		set_intensity(out, h, k, l, val1/val2);
		add_item(icommon, h, k, l);
414
415

	}
416
	ncom = num_items(icommon);
Thomas White's avatar
Thomas White committed
417
418
	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);
419

Thomas White's avatar
Thomas White committed
420
421
	STATUS("%i,%i reflections: %i in common\n",
	       num_items(i1), num_items(i2), ncom);
Thomas White's avatar
Thomas White committed
422

423
	R1 = stat_r1_ignore(ref1, ref2_transformed, icommon, &scale_r1fi);
424
	STATUS("R1(F) = %5.4f %% (scale=%5.2e) (ignoring negative intensities)\n",
425
	       R1*100.0, scale_r1fi);
Thomas White's avatar
Thomas White committed
426

427
	R1 = stat_r1_zero(ref1, ref2_transformed, icommon, &scale_r1);
428
	STATUS("R1(F) = %5.4f %% (scale=%5.2e) (zeroing negative intensities)\n",
429
	       R1*100.0, scale_r1);
Thomas White's avatar
Thomas White committed
430

Thomas White's avatar
Thomas White committed
431
	R2 = stat_r2(ref1, ref2_transformed, icommon, &scale_r2);
432
	STATUS("R2(I) = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale_r2);
Thomas White's avatar
Thomas White committed
433

434
	R1i = stat_r1_i(ref1, ref2_transformed, icommon, &scale_r1i);
435
	STATUS("R1(I) = %5.4f %% (scale=%5.2e)\n", R1i*100.0, scale_r1i);
Thomas White's avatar
Thomas White committed
436

437
	Rdiff = stat_rdiff_ignore(ref1, ref2_transformed, icommon, &scale_rdig);
Thomas White's avatar
Thomas White committed
438
	STATUS("Rint(F) = %5.4f %% (scale=%5.2e) (ignoring negative intensities)\n",
439
	       Rdiff*100.0, scale_rdig);
Thomas White's avatar
Thomas White committed
440
441

	Rdiff = stat_rdiff_zero(ref1, ref2_transformed, icommon, &scale);
Thomas White's avatar
Thomas White committed
442
443
444
445
446
447
	STATUS("Rint(F) = %5.4f %% (scale=%5.2e) (zeroing negative intensities)\n",
	       Rdiff*100.0, scale);

	Rdiff = stat_rdiff_intensity(ref1, ref2_transformed, icommon,
	                             &scale_rintint);
	STATUS("Rint(I) = %5.4f %% (scale=%5.2e)\n",
448
	       Rdiff*100.0, scale_rintint);
Thomas White's avatar
Thomas White committed
449
450
451
452
453
454
455
456
457
458
459

	pearson = stat_pearson_i(ref1, ref2_transformed, icommon);
	STATUS("Pearson r(I) = %5.4f\n", pearson);

	pearson = stat_pearson_f_ignore(ref1, ref2_transformed, icommon);
	STATUS("Pearson r(F) = %5.4f (ignoring negative intensities)\n",
	       pearson);

	pearson = stat_pearson_f_zero(ref1, ref2_transformed, icommon);
	STATUS("Pearson r(F) = %5.4f (zeroing negative intensities)\n",
	       pearson);
Thomas White's avatar
Thomas White committed
460

Thomas White's avatar
Thomas White committed
461
	if ( config_shells ) {
462
		plot_shells(ref1, ref2_transformed, icommon, scale_r1fi,
463
		            cell, sym, rmin_fix, rmax_fix);
Thomas White's avatar
Thomas White committed
464
465
	}

466
	if ( outfile != NULL ) {
467

468
469
		write_reflections(outfile, icommon, out, NULL,
		                  NULL, NULL, cell);
470
471
		STATUS("Sigma(I) values in output file are not meaningful.\n");

472
	}
473
474
475

	return 0;
}