compare_hkl.c 10.3 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
26
/*
 * 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 "sfac.h"
#include "reflections.h"
Thomas White's avatar
Thomas White committed
27
#include "statistics.h"
28
#include "symmetry.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
42
	printf(
"Compare intensity lists.\n"
"\n"
"  -h, --help                 Display this help message.\n"
"  -o, --output=<filename>    Specify output filename for correction factor.\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
static void plot_shells(const double *ref1, const double *ref2,
53
                        ReflItemList *items, double scale, UnitCell *cell,
54
                        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

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

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

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

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

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

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

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

	}
110

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

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

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

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

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

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

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

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

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

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

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

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

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

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

	}

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

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

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

	}

	fclose(fh);
}


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

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

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

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

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

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

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

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

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

	}

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

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

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

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

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

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

Thomas White's avatar
Thomas White committed
355
356
	/* List for output scale factor map */
	out = new_list_intensity();
357

358
359
360
361
362
	/* 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
363
364
		struct refl_item *it;
		signed int h, k, l;
365
366
		signed int he, ke, le;
		double val1, val2;
367
368
		double sig1, sig2;
		int ig = 0;
369
		double d;
370

371
		it = get_item(i1, i);
Thomas White's avatar
Thomas White committed
372
		h = it->h;  k = it->k;  l = it->l;
373

374
		if ( !find_unique_equiv(i2, h, k, l, sym, &he, &ke, &le) ) {
Thomas White's avatar
Thomas White committed
375
376
			//STATUS("%i %i %i not matched (%f nm).\n", h, k, l,
			//       1.0/(2.0*resolution(cell, h, k, l)/1e9));
377
378
			continue;
		}
379

380
381
		val1 = lookup_intensity(ref1, h, k, l);
		val2 = lookup_intensity(ref2, he, ke, le);
382
		sig1 = lookup_sigma(esd1, h, k, l);
383
		sig2 = lookup_sigma(esd2, he, ke, le);
384
385
386
387
388
389
390
391
392

		if ( val1 < 3.0 * sig1 ) {
			rej1++;
			ig = 1;
		}
		if ( val2 < 3.0 * sig2 ) {
			rej2++;
			ig = 1;
		}
393
394
395
396
397

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

398
		//if ( ig ) continue;
399

400
401
402
		set_intensity(ref2_transformed, h, k, l, val2);
		set_intensity(out, h, k, l, val1/val2);
		add_item(icommon, h, k, l);
403
404

	}
405
	ncom = num_items(icommon);
406
407
408
409
	STATUS("%i reflections with I < 3.0*sigma(I) rejected from '%s'\n",
	       rej1, afile);
	STATUS("%i reflections with I < 3.0*sigma(I) rejected from '%s'\n",
	       rej2, bfile);
410

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

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

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

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

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

428
	Rdiff = stat_rdiff_ignore(ref1, ref2_transformed, icommon, &scale_rdig);
Thomas White's avatar
Thomas White committed
429
	STATUS("Rint(F) = %5.4f %% (scale=%5.2e) (ignoring negative intensities)\n",
430
	       Rdiff*100.0, scale_rdig);
Thomas White's avatar
Thomas White committed
431
432

	Rdiff = stat_rdiff_zero(ref1, ref2_transformed, icommon, &scale);
Thomas White's avatar
Thomas White committed
433
434
435
436
437
438
	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",
439
	       Rdiff*100.0, scale_rintint);
Thomas White's avatar
Thomas White committed
440
441
442
443
444
445
446
447
448
449
450

	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
451

Thomas White's avatar
Thomas White committed
452
	if ( config_shells ) {
453
		plot_shells(ref1, ref2_transformed, icommon, scale_r1fi,
454
		            cell, sym, rmin_fix, rmax_fix);
Thomas White's avatar
Thomas White committed
455
456
	}

457
	if ( outfile != NULL ) {
458

459
460
		write_reflections(outfile, icommon, out, NULL,
		                  NULL, NULL, cell);
461
462
		STATUS("Sigma(I) values in output file are not meaningful.\n");

463
	}
464
465
466

	return 0;
}