compare_hkl.c 11.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
                        ReflItemList *items, double scale, UnitCell *cell,
                        const char *sym, ReflItemList *characterise,
53
                        unsigned int *char_counts, const double *sigma)
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
72
	double snr_total = 0;
	int nmeas = 0;
73
	int nmeastot = 0;
Thomas White's avatar
Thomas White committed
74
75

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

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

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

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

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

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

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

	}
111
112
113
114
115
116

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

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

117
118
119
120
	/* FIXME: Fixed resolution shells */
	rmin = 0.120e9;
	rmax = 1.172e9;

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
		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 */
Thomas White's avatar
Thomas White committed
176
		if ( (h==0) && (k==0) && (l%2) ) continue;
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

		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 */
Thomas White's avatar
Thomas White committed
213
		if ( (h==0) && (k==0) && (l%2) ) continue;
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230

		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);
231
232
233
234
235
		snr[bin] += (lookup_intensity(ref1, h, k, l) /
		                              lookup_intensity(sigma, h, k, l));
		snr_total += (lookup_intensity(ref1, h, k, l) /
		                              lookup_intensity(sigma, h, k, l));
		nmeas++;
236
		nmeastot += lookup_count(char_counts, h, k, l);
237
238

	}
239
	STATUS("overall <snr> = %f\n", snr_total/(double)nmeas);
240
241
	STATUS("%i measurements in total.\n", nmeastot);
	STATUS("%i reflections in total.\n", nmeas);
Thomas White's avatar
Thomas White committed
242

Thomas White's avatar
Thomas White committed
243
	den = 0.0;
Thomas White's avatar
Thomas White committed
244
	ctot = 0;
Thomas White's avatar
Thomas White committed
245
246
247
248
	for ( i=0; i<num_items(items); i++ ) {

		struct refl_item *it;
		signed int h, k, l;
249
		double d;
Thomas White's avatar
Thomas White committed
250
		int bin;
251
		double i1, i2, f1, f2;
252
		int j;
Thomas White's avatar
Thomas White committed
253
254
255
256

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

257
		/* FIXME: Reflection condition */
Thomas White's avatar
Thomas White committed
258
		if ( (h==0) && (k==0) && (l%2) ) continue;
259

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

262
263
264
265
266
267
268
269
270
271
272
273
		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
274
275

		i1 = lookup_intensity(ref1, h, k, l);
Thomas White's avatar
Thomas White committed
276
277
		//if ( i1 < 0.0 ) continue;
		//f1 = sqrt(i1);
278
		i2 = lookup_intensity(ref2, h, k, l);
Thomas White's avatar
Thomas White committed
279
280
281
		//if ( i2 < 0.0 ) continue;
		//f2 = sqrt(i2);
		i2 *= scale;
Thomas White's avatar
Thomas White committed
282

Thomas White's avatar
Thomas White committed
283
284
		num[bin] += fabs(i1 - i2);
		den += i1;// + i2) / 2.0;
Thomas White's avatar
Thomas White committed
285
		ctot++;
286
		cts[bin]++;
Thomas White's avatar
Thomas White committed
287
288
289
290
291

	}

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

292
		double r, cen;
293
		cen = rmins[i] + (rmaxs[i] - rmins[i])/2.0;
294
		r = (num[i]/den)*((double)ctot/cts[i]);
295
		fprintf(fh, "%f %f %i %i %i %f %f\n", cen*1.0e-9, r*100.0,
296
297
		                            measured[i],
		                            possible[i], measurements[i],
298
299
		                            (float)measurements[i]/measured[i],
		                            (snr[i]/(double)measured[i]));
Thomas White's avatar
Thomas White committed
300
301
302
303
304
305
306

	}

	fclose(fh);
}


307
308
309
310
311
int main(int argc, char *argv[])
{
	int c;
	double *ref1;
	double *ref2;
312
	double *ref2_transformed;
313
	double *out;
314
	UnitCell *cell;
315
316
317
	char *outfile = NULL;
	char *afile = NULL;
	char *bfile = NULL;
Thomas White's avatar
Thomas White committed
318
	char *sym = NULL;
319
	double scale, scale_r2, scale_rdig, R1, R2, R1i, Rdiff, pearson;
320
	double scale_rintint, scale_r1i, scale_r1, scale_r1fi;
Thomas White's avatar
Thomas White committed
321
322
	int i, ncom;
	ReflItemList *i1, *i2, *icommon;
Thomas White's avatar
Thomas White committed
323
	int config_shells = 0;
324
	char *pdb = NULL;
325
326
327
328
	double *esd1;
	double *esd2;
	int rej1 = 0;
	int rej2 = 0;
329
	unsigned int *cts1;
330
331
332
333
334

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"output",             1, NULL,               'o'},
335
		{"symmetry",           1, NULL,               'y'},
Thomas White's avatar
Thomas White committed
336
		{"shells",             0, &config_shells,     1},
337
		{"pdb",                1, NULL,               'p'},
338
339
340
341
		{0, 0, NULL, 0}
	};

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

		switch (c) {
Thomas White's avatar
Thomas White committed
345
		case 'h' :
346
347
348
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
349
		case 'o' :
350
351
352
			outfile = strdup(optarg);
			break;

353
354
355
356
		case 'y' :
			sym = strdup(optarg);
			break;

357
358
359
360
		case 'p' :
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
361
		case 0 :
362
363
			break;

Thomas White's avatar
Thomas White committed
364
		default :
365
366
367
368
369
			return 1;
		}

	}

370
371
372
373
374
	if ( argc != (optind+2) ) {
		ERROR("Please provide exactly two HKL files to compare.\n");
		return 1;
	}

375
376
377
378
	if ( sym == NULL ) {
		sym = strdup("1");
	}

379
380
381
	afile = strdup(argv[optind++]);
	bfile = strdup(argv[optind]);

382
383
384
385
386
387
388
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

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

Thomas White's avatar
Thomas White committed
389
	ref1 = new_list_intensity();
390
	esd1 = new_list_sigma();
391
392
	cts1 = new_list_count();
	i1 = read_reflections(afile, ref1, NULL, cts1, esd1);
Thomas White's avatar
Thomas White committed
393
	if ( i1 == NULL ) {
394
395
396
		ERROR("Couldn't open file '%s'\n", afile);
		return 1;
	}
Thomas White's avatar
Thomas White committed
397
	ref2 = new_list_intensity();
398
399
	esd2 = new_list_sigma();
	i2 = read_reflections(bfile, ref2, NULL, NULL, esd2);
Thomas White's avatar
Thomas White committed
400
	if ( i2 == NULL ) {
401
402
403
404
		ERROR("Couldn't open file '%s'\n", bfile);
		return 1;
	}

Thomas White's avatar
Thomas White committed
405
406
	/* List for output scale factor map */
	out = new_list_intensity();
407

408
409
410
411
412
	/* 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
413
414
		struct refl_item *it;
		signed int h, k, l;
415
416
		signed int he, ke, le;
		double val1, val2;
417
418
		double sig1, sig2;
		int ig = 0;
419
		double d;
420

421
		it = get_item(i1, i);
Thomas White's avatar
Thomas White committed
422
		h = it->h;  k = it->k;  l = it->l;
423

424
		if ( !find_unique_equiv(i2, h, k, l, sym, &he, &ke, &le) ) {
Thomas White's avatar
Thomas White committed
425
426
			//STATUS("%i %i %i not matched (%f nm).\n", h, k, l,
			//       1.0/(2.0*resolution(cell, h, k, l)/1e9));
427
428
			continue;
		}
429

430
431
		val1 = lookup_intensity(ref1, h, k, l);
		val2 = lookup_intensity(ref2, he, ke, le);
432
		sig1 = lookup_sigma(esd1, h, k, l);
433
		sig2 = lookup_sigma(esd2, he, ke, le);
434
435
436
437
438
439
440
441
442

		if ( val1 < 3.0 * sig1 ) {
			rej1++;
			ig = 1;
		}
		if ( val2 < 3.0 * sig2 ) {
			rej2++;
			ig = 1;
		}
443
444
445
446
447

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

448
		//if ( ig ) continue;
449

450
451
452
		set_intensity(ref2_transformed, h, k, l, val2);
		set_intensity(out, h, k, l, val1/val2);
		add_item(icommon, h, k, l);
453
454

	}
455
	ncom = num_items(icommon);
456
457
458
459
	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);
460

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

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

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

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

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

478
	Rdiff = stat_rdiff_ignore(ref1, ref2_transformed, icommon, &scale_rdig);
Thomas White's avatar
Thomas White committed
479
	STATUS("Rint(F) = %5.4f %% (scale=%5.2e) (ignoring negative intensities)\n",
480
	       Rdiff*100.0, scale_rdig);
Thomas White's avatar
Thomas White committed
481
482

	Rdiff = stat_rdiff_zero(ref1, ref2_transformed, icommon, &scale);
Thomas White's avatar
Thomas White committed
483
484
485
486
487
488
	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",
489
	       Rdiff*100.0, scale_rintint);
Thomas White's avatar
Thomas White committed
490
491
492
493
494
495
496
497
498
499
500

	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
501

Thomas White's avatar
Thomas White committed
502
	if ( config_shells ) {
503
504
		plot_shells(ref1, ref2_transformed, icommon, scale_r1fi,
		            cell, sym, i1, cts1, esd1);
Thomas White's avatar
Thomas White committed
505
506
	}

507
	if ( outfile != NULL ) {
508
509
510
511

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

512
	}
513
514
515

	return 0;
}