compare_hkl.c 10.8 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
48
49
"\n");
}


Thomas White's avatar
Thomas White committed
50
static void plot_shells(const double *ref1, const double *ref2,
51
52
53
                        ReflItemList *items, double scale, UnitCell *cell,
                        const char *sym, ReflItemList *characterise,
                        unsigned int *char_counts)
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];
Thomas White's avatar
Thomas White committed
64
	double den;
65
66
	double rmin, rmax;
	signed int h, k, l;
Thomas White's avatar
Thomas White committed
67
	int i;
Thomas White's avatar
Thomas White committed
68
	int ctot;
Thomas White's avatar
Thomas White committed
69
70
71
	FILE *fh;

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

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

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

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

		struct refl_item *it;
		signed int h, k, l;
96
		double d;
Thomas White's avatar
Thomas White committed
97
98
99
100

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

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

	}
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
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
161
162
163
164
165
166
167
168
169
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
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218

	STATUS("%f -> %f\n", rmin/1e9, rmax/1e9);

	/* Increase the max just a little bit */
	rmax += 0.001e9;

	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);

		rmaxs[i-1] = r;
		rmins[i] = r;

		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);

#if 0
	/* FIXME: Fixed resolution shells */
	rmins[0] = 0.121065;  rmaxs[0] = 0.552486;
	rmins[1] = 0.552486;  rmaxs[1] = 0.690186;
	rmins[2] = 0.690186;  rmaxs[2] = 0.787787;
	rmins[3] = 0.787787;  rmaxs[3] = 0.865813;
	rmins[4] = 0.865813;  rmaxs[4] = 0.931853;
	rmins[5] = 0.931853;  rmaxs[5] = 0.989663;
	rmins[6] = 0.989663;  rmaxs[6] = 1.041409;
	rmins[7] = 1.041409;  rmaxs[7] = 1.088467;
	rmins[8] = 1.088467;  rmaxs[8] = 1.131775;
	rmins[9] = 1.131775;  rmaxs[9] = 1.172000;
	for ( i=0; i<NBINS; i++ ) {
		rmins[i] *= 1e9;
		rmaxs[i] *= 1e9;
	}
#endif

	/* 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;

		/* FIXME: Reflection condition */
		// if ( (h==0) && (k==0) && (l%2) ) continue;

		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);

	/* Characterise the first data set (only) */
	for ( i=0; i<num_items(characterise); i++ ) {

		struct refl_item *it;
		signed int h, k, l;
		double d;
		int bin;
		int j;

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

		/* FIXME: Reflection condition */
		//if ( (h==0) && (k==0) && (l%2) ) continue;

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

		bin = -1;
		for ( j=0; j<NBINS; j++ ) {
			if ( (d>rmins[j]) && (d<=rmaxs[j]) ) {
				bin = j;
				break;
			}
		}
		if ( bin == -1 ) {
			ERROR("Warnung! %i %i %i %f\n", h, k, l, d/1e9);
			continue;
		}

		measured[bin]++;
		measurements[bin] += lookup_count(char_counts, h, k, l);

	}
Thomas White's avatar
Thomas White committed
219

Thomas White's avatar
Thomas White committed
220
	den = 0.0;
Thomas White's avatar
Thomas White committed
221
	ctot = 0;
Thomas White's avatar
Thomas White committed
222
223
224
225
	for ( i=0; i<num_items(items); i++ ) {

		struct refl_item *it;
		signed int h, k, l;
226
		double d;
Thomas White's avatar
Thomas White committed
227
		int bin;
228
		double i1, i2, f1, f2;
229
		int j;
Thomas White's avatar
Thomas White committed
230
231
232
233

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

234
235
236
		/* FIXME: Reflection condition */
		//if ( (h==0) && (k==0) && (l%2) ) continue;

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

239
240
241
242
243
244
245
246
247
248
249
250
		bin = -1;
		for ( j=0; j<NBINS; j++ ) {
			if ( (d>rmins[j]) && (d<=rmaxs[j]) ) {
				bin = j;
				break;
			}
		}
		if ( bin == -1 ) {
			ERROR("Warnung! %i %i %i %f\n", h, k, l, d/1e9);
			abort();
			continue;
		}
Thomas White's avatar
Thomas White committed
251
252

		i1 = lookup_intensity(ref1, h, k, l);
253
254
255
256
257
		if ( i1 < 0.0 ) continue;
		f1 = sqrt(i1);
		i2 = lookup_intensity(ref2, h, k, l);
		if ( i2 < 0.0 ) continue;
		f2 = sqrt(i2);
258
		f2 *= scale;
Thomas White's avatar
Thomas White committed
259

260
		num[bin] += fabs(f1 - f2);
Thomas White's avatar
Thomas White committed
261
		den += (f1 + f2) / 2.0;
Thomas White's avatar
Thomas White committed
262
		ctot++;
263
		cts[bin]++;
Thomas White's avatar
Thomas White committed
264
265
266
267
268

	}

	for ( i=0; i<NBINS; i++ ) {

269
		double r, cen;
270
		cen = rmins[i] + (rmaxs[i] - rmins[i])/2.0;
271
		r = (num[i]/den)*((double)ctot/cts[i]);
272
273
274
275
		fprintf(fh, "%f %f %i %i %i %f\n", cen*1.0e-9, r*100.0,
		                            measured[i],
		                            possible[i], measurements[i],
		                            (float)measurements[i]/measured[i]);
Thomas White's avatar
Thomas White committed
276
277
278
279
280
281
282

	}

	fclose(fh);
}


283
284
285
286
287
int main(int argc, char *argv[])
{
	int c;
	double *ref1;
	double *ref2;
288
	double *ref2_transformed;
289
	double *out;
290
	UnitCell *cell;
291
292
293
	char *outfile = NULL;
	char *afile = NULL;
	char *bfile = NULL;
Thomas White's avatar
Thomas White committed
294
	char *sym = NULL;
295
	double scale, scale_r2, scale_rdig, R1, R2, R1i, Rdiff, pearson;
Thomas White's avatar
Thomas White committed
296
	double scale_rintint;
Thomas White's avatar
Thomas White committed
297
298
	int i, ncom;
	ReflItemList *i1, *i2, *icommon;
Thomas White's avatar
Thomas White committed
299
	int config_shells = 0;
300
	char *pdb = NULL;
301
302
303
304
	double *esd1;
	double *esd2;
	int rej1 = 0;
	int rej2 = 0;
305
	unsigned int *cts1;
306
307
308
309
310

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"output",             1, NULL,               'o'},
311
		{"symmetry",           1, NULL,               'y'},
Thomas White's avatar
Thomas White committed
312
		{"shells",             0, &config_shells,     1},
313
		{"pdb",                1, NULL,               'p'},
314
315
316
317
		{0, 0, NULL, 0}
	};

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

		switch (c) {
Thomas White's avatar
Thomas White committed
321
		case 'h' :
322
323
324
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
325
		case 'o' :
326
327
328
			outfile = strdup(optarg);
			break;

329
330
331
332
		case 'y' :
			sym = strdup(optarg);
			break;

333
334
335
336
		case 'p' :
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
337
		case 0 :
338
339
			break;

Thomas White's avatar
Thomas White committed
340
		default :
341
342
343
344
345
			return 1;
		}

	}

346
347
348
349
350
	if ( argc != (optind+2) ) {
		ERROR("Please provide exactly two HKL files to compare.\n");
		return 1;
	}

351
352
353
354
	if ( sym == NULL ) {
		sym = strdup("1");
	}

355
356
357
	afile = strdup(argv[optind++]);
	bfile = strdup(argv[optind]);

358
359
360
361
362
363
364
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

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

Thomas White's avatar
Thomas White committed
365
	ref1 = new_list_intensity();
366
	esd1 = new_list_sigma();
367
368
	cts1 = new_list_count();
	i1 = read_reflections(afile, ref1, NULL, cts1, esd1);
Thomas White's avatar
Thomas White committed
369
	if ( i1 == NULL ) {
370
371
372
		ERROR("Couldn't open file '%s'\n", afile);
		return 1;
	}
Thomas White's avatar
Thomas White committed
373
	ref2 = new_list_intensity();
374
375
	esd2 = new_list_sigma();
	i2 = read_reflections(bfile, ref2, NULL, NULL, esd2);
Thomas White's avatar
Thomas White committed
376
	if ( i2 == NULL ) {
377
378
379
380
		ERROR("Couldn't open file '%s'\n", bfile);
		return 1;
	}

Thomas White's avatar
Thomas White committed
381
382
	/* List for output scale factor map */
	out = new_list_intensity();
383

384
385
386
387
388
	/* 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
389
390
		struct refl_item *it;
		signed int h, k, l;
391
392
		signed int he, ke, le;
		double val1, val2;
393
394
		double sig1, sig2;
		int ig = 0;
395
		double d;
396

397
		it = get_item(i1, i);
Thomas White's avatar
Thomas White committed
398
		h = it->h;  k = it->k;  l = it->l;
399

400
		if ( !find_unique_equiv(i2, h, k, l, sym, &he, &ke, &le) ) {
Thomas White's avatar
Thomas White committed
401
402
			//STATUS("%i %i %i not matched (%f nm).\n", h, k, l,
			//       1.0/(2.0*resolution(cell, h, k, l)/1e9));
403
404
			continue;
		}
405

406
407
		val1 = lookup_intensity(ref1, h, k, l);
		val2 = lookup_intensity(ref2, he, ke, le);
408
		sig1 = lookup_sigma(esd1, h, k, l);
409
		sig2 = lookup_sigma(esd2, he, ke, le);
410
411
412
413
414
415
416
417
418

		if ( val1 < 3.0 * sig1 ) {
			rej1++;
			ig = 1;
		}
		if ( val2 < 3.0 * sig2 ) {
			rej2++;
			ig = 1;
		}
419
420
421
422
423

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

424
		//if ( ig ) continue;
425

426
427
428
		set_intensity(ref2_transformed, h, k, l, val2);
		set_intensity(out, h, k, l, val1/val2);
		add_item(icommon, h, k, l);
429
430

	}
431
	ncom = num_items(icommon);
432
433
434
435
	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);
436

Thomas White's avatar
Thomas White committed
437
438
	STATUS("%i,%i reflections: %i in common\n",
	       num_items(i1), num_items(i2), ncom);
Thomas White's avatar
Thomas White committed
439
440

	R1 = stat_r1_ignore(ref1, ref2_transformed, icommon, &scale);
441
	STATUS("R1(F) = %5.4f %% (scale=%5.2e) (ignoring negative intensities)\n",
Thomas White's avatar
Thomas White committed
442
443
444
	       R1*100.0, scale);

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

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

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

454
	Rdiff = stat_rdiff_ignore(ref1, ref2_transformed, icommon, &scale_rdig);
Thomas White's avatar
Thomas White committed
455
	STATUS("Rint(F) = %5.4f %% (scale=%5.2e) (ignoring negative intensities)\n",
Thomas White's avatar
Thomas White committed
456
457
458
	       Rdiff*100.0, scale);

	Rdiff = stat_rdiff_zero(ref1, ref2_transformed, icommon, &scale);
Thomas White's avatar
Thomas White committed
459
460
461
462
463
464
	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",
Thomas White's avatar
Thomas White committed
465
466
467
468
469
470
471
472
473
474
475
476
	       Rdiff*100.0, scale);

	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
477

Thomas White's avatar
Thomas White committed
478
	if ( config_shells ) {
479
480
		plot_shells(ref1, ref2_transformed, icommon, scale_rdig, cell,
		            sym, i1, cts1);
Thomas White's avatar
Thomas White committed
481
482
	}

483
	if ( outfile != NULL ) {
Thomas White's avatar
Thomas White committed
484
		write_reflections(outfile, icommon, out, NULL, NULL, cell);
485
	}
486
487
488

	return 0;
}