compare_hkl.c 11 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
#include <gsl/gsl_errno.h>
25
26

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


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


36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
enum r_shell
{
	R_SHELL_NONE,
	R_SHELL_R1I,
	R_SHELL_R1F,
	R_SHELL_RSPLIT,
};


static enum r_shell get_r_shell(const char *s)
{
	if ( strcmp(s, "r1i") == 0 ) return R_SHELL_R1I;
	if ( strcmp(s, "r1f") == 0 ) return R_SHELL_R1F;
	if ( strcmp(s, "rsplit") == 0 ) return R_SHELL_RSPLIT;

	ERROR("Unknown R-factor '%s' - try '--shells=rsplit', or --help for"
	      " more possibilities.\n", s);
	exit(1);
}


57
58
static void show_help(const char *s)
{
59
	printf("Syntax: %s [options] <file1.hkl> <file2.hkl>\n\n", s);
60
61
62
63
	printf(
"Compare intensity lists.\n"
"\n"
"  -h, --help                 Display this help message.\n"
Thomas White's avatar
Thomas White committed
64
"  -o, --ratio=<filename>     Specify output filename for ratios.\n"
65
"  -y, --symmetry=<sym>       The symmetry of both the input files.\n"
66
"  -p, --pdb=<filename>       PDB file to use (default: molecule.pdb).\n"
67
68
"      --shells=<FoM>         Plot this figure of merit in resolution shells.\n"
"                              Choose from: 'Rsplit', 'R1f' and 'R1i'.\n"
69
70
"      --rmin=<res>           Fix lower resolution limit for --shells (m^-1).\n"
"      --rmax=<res>           Fix upper resolution limit for --shells (m^-1).\n"
71
72
73
74
"\n");
}


75
static void plot_shells(RefList *list1, RefList *list2, double scale,
76
77
                        UnitCell *cell, double rmin_fix, double rmax_fix,
                        enum r_shell config_shells)
Thomas White's avatar
Thomas White committed
78
79
{
	double num[NBINS];
Thomas White's avatar
Thomas White committed
80
	int cts[NBINS];
81
82
83
84
85
	unsigned int measurements[NBINS];
	unsigned int measured[NBINS];
	double total_vol, vol_per_shell;
	double rmins[NBINS];
	double rmaxs[NBINS];
86
	double snr[NBINS];
87
	double rmin, rmax;
Thomas White's avatar
Thomas White committed
88
	int i;
Thomas White's avatar
Thomas White committed
89
90
	Reflection *refl1;
	RefListIterator *iter;
91
92
93
	FILE *fh;
	double den;
	int ctot, nout;
Thomas White's avatar
Thomas White committed
94
95

	if ( cell == NULL ) {
Thomas White's avatar
Thomas White committed
96
		ERROR("Need the unit cell to plot resolution shells.\n");
Thomas White's avatar
Thomas White committed
97
98
99
100
101
		return;
	}

	for ( i=0; i<NBINS; i++ ) {
		num[i] = 0.0;
Thomas White's avatar
Thomas White committed
102
		cts[i] = 0;
103
104
		measured[i] = 0;
		measurements[i] = 0;
105
		snr[i] = 0;
Thomas White's avatar
Thomas White committed
106
107
	}

108
109
	/* Find resolution limits */
	resolution_limits(list1, cell, &rmin, &rmax);
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

Thomas White's avatar
Thomas White committed
120
	/* Calculate the resolution bins */
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
		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);

147
	den = 0.0;  ctot = 0;  nout = 0;
Thomas White's avatar
Thomas White committed
148
149
	for ( refl1 = first_refl(list1, &iter);
	      refl1 != NULL;
150
151
	      refl1 = next_refl(refl1, iter) )
	{
Thomas White's avatar
Thomas White committed
152
		signed int h, k, l;
153
		double d;
Thomas White's avatar
Thomas White committed
154
		int bin;
155
		double i1, i2, f1, f2;
156
		int j;
157
		Reflection *refl2;
Thomas White's avatar
Thomas White committed
158

Thomas White's avatar
Thomas White committed
159
		get_indices(refl1, &h, &k, &l);
160
		d = 2.0 * resolution(cell, h, k, l);
Thomas White's avatar
Thomas White committed
161

162
163
164
165
166
167
168
		bin = -1;
		for ( j=0; j<NBINS; j++ ) {
			if ( (d>rmins[j]) && (d<=rmaxs[j]) ) {
				bin = j;
				break;
			}
		}
169
170

		/* Outside resolution range? */
171
172
173
174
		if ( bin == -1 ) {
			nout++;
			continue;
		}
Thomas White's avatar
Thomas White committed
175

176
177
178
		refl2 = find_refl(list2, h, k, l);
		if ( refl2 == NULL ) continue;

Thomas White's avatar
Thomas White committed
179
		i1 = get_intensity(refl1);
180
181
182
183
184
185
186
187
		i2 = get_intensity(refl2);
		f1 = i1 > 0.0 ? sqrt(i1) : 0.0;
		f2 = i2 > 0.0 ? sqrt(i2) : 0.0;

		switch ( config_shells ) {

			case R_SHELL_RSPLIT :
			num[bin] += fabs(i1 - scale*i2);
188
			den += i1 + scale*i2;
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
			break;

			case R_SHELL_R1I :
			num[bin] += fabs(i1 - scale*i2);
			den += i1;
			break;

			case R_SHELL_R1F :
			num[bin] += fabs(f1 - scale*f2);
			den += f1;
			break;

			default : break;

		}
Thomas White's avatar
Thomas White committed
204

Thomas White's avatar
Thomas White committed
205
		ctot++;
206
		cts[bin]++;
Thomas White's avatar
Thomas White committed
207
208
	}

209
210
211
212
213
	if ( nout ) {
		STATUS("Warning; %i reflections outside resolution range.\n",
		       nout);
	}

214
215
216
217
218
219
	fh = fopen("shells.dat", "w");
	if ( fh == NULL ) {
		ERROR("Couldn't open 'shells.dat'\n");
		return;
	}

220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
	switch ( config_shells ) {

		case R_SHELL_RSPLIT :
		fprintf(fh, "1/d centre   Rsplit / %%\n");
		break;

		case R_SHELL_R1I :
		fprintf(fh, "1/d centre   R1(I) / %%\n");
		break;

		case R_SHELL_R1F :
		fprintf(fh, "1/d centre   R1(F) ignoring -ves / %%\n");
		break;

		default :
		fprintf(fh, "1/d centre   0.0\n");
		break;

	}
239

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

242
		double r, cen;
243
		cen = rmins[i] + (rmaxs[i] - rmins[i])/2.0;
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261

		switch ( config_shells ) {

			case R_SHELL_RSPLIT :
			r = (2.0*(num[i]/den)*((double)ctot/cts[i]))/sqrt(2.0);
			break;

			case R_SHELL_R1I :
			case R_SHELL_R1F :
			r = (num[i]/den) * ((double)ctot/cts[i]);
			break;

			default :
			r = 0.0;
			break;

		}

262
263
		fprintf(fh, "%10.3f %10.2f %10i\n",
		        cen*1.0e-9, r*100.0, cts[i]);
Thomas White's avatar
Thomas White committed
264
265
266
267

	}

	fclose(fh);
Thomas White's avatar
Thomas White committed
268
269

	STATUS("Resolution shell information written to shells.dat.\n");
Thomas White's avatar
Thomas White committed
270
271
272
}


273
274
275
int main(int argc, char *argv[])
{
	int c;
276
	UnitCell *cell;
Thomas White's avatar
Thomas White committed
277
	char *ratiofile = NULL;
278
279
	char *afile = NULL;
	char *bfile = NULL;
Thomas White's avatar
Thomas White committed
280
281
	char *sym_str = NULL;
	SymOpList *sym;
282
	double scale, scale_r2, scale_rdig, R1, R2, R1i, Rdiff, pearson;
283
	double scale_rintint, scale_r1i, scale_r1, scale_r1fi;
Thomas White's avatar
Thomas White committed
284
285
286
	int ncom;
	RefList *list1;
	RefList *list2;
287
288
	RefList *list1_raw;
	RefList *list2_raw;
Thomas White's avatar
Thomas White committed
289
	RefList *ratio;
290
	enum r_shell config_shells = R_SHELL_NONE;
291
	char *pdb = NULL;
292
293
	float rmin_fix = -1.0;
	float rmax_fix = -1.0;
Thomas White's avatar
Thomas White committed
294
295
	Reflection *refl1;
	RefListIterator *iter;
296
	int config_unity = 0;
297
	double scale_for_shells;
298
299
300
301

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
Thomas White's avatar
Thomas White committed
302
		{"ratio" ,             1, NULL,               'o'},
303
		{"symmetry",           1, NULL,               'y'},
304
		{"shells",             1, NULL,               4},
305
		{"pdb",                1, NULL,               'p'},
306
307
		{"rmin",               1, NULL,               2},
		{"rmax",               1, NULL,               3},
308
309
310
311
		{0, 0, NULL, 0}
	};

	/* Short options */
312
313
314
	while ((c = getopt_long(argc, argv, "ho:y:p:u",
	                        longopts, NULL)) != -1)
	{
315
316

		switch (c) {
Thomas White's avatar
Thomas White committed
317
		case 'h' :
318
319
320
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
321
		case 'o' :
Thomas White's avatar
Thomas White committed
322
			ratiofile = strdup(optarg);
323
324
			break;

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

329
330
331
332
		case 'p' :
			pdb = strdup(optarg);
			break;

333
334
335
336
		case 'u' :
			config_unity = 1;
			break;

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

340
341
342
343
344
345
346
347
348
349
350
351
352
353
		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;

354
355
356
357
		case 4 :
			config_shells = get_r_shell(optarg);
			break;

Thomas White's avatar
Thomas White committed
358
		default :
359
360
361
362
363
			return 1;
		}

	}

364
365
366
367
368
	if ( argc != (optind+2) ) {
		ERROR("Please provide exactly two HKL files to compare.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
369
370
	if ( sym_str == NULL ) {
		sym_str = strdup("1");
371
	}
Thomas White's avatar
Thomas White committed
372
373
	sym = get_pointgroup(sym_str);
	free(sym_str);
374

375
376
377
	afile = strdup(argv[optind++]);
	bfile = strdup(argv[optind]);

378
379
380
381
382
383
384
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

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

385
386
	list1_raw = read_reflections(afile);
	if ( list1_raw == NULL ) {
Thomas White's avatar
Thomas White committed
387
		ERROR("Couldn't read file '%s'\n", afile);
388
389
		return 1;
	}
Thomas White's avatar
Thomas White committed
390

391
392
	list2_raw = read_reflections(bfile);
	if ( list2_raw == NULL ) {
Thomas White's avatar
Thomas White committed
393
		ERROR("Couldn't read file '%s'\n", bfile);
394
395
396
		return 1;
	}

Thomas White's avatar
Thomas White committed
397
	/* Check that the intensities have the correct symmetry */
398
	if ( check_list_symmetry(list1_raw, sym) ) {
Thomas White's avatar
Thomas White committed
399
		ERROR("The first input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
400
		      " have symmetry %s\n", symmetry_name(sym));
Thomas White's avatar
Thomas White committed
401
402
		return 1;
	}
403
	if ( check_list_symmetry(list2_raw, sym) ) {
Thomas White's avatar
Thomas White committed
404
		ERROR("The second input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
405
		      " have symmetry %s\n", symmetry_name(sym));
Thomas White's avatar
Thomas White committed
406
407
408
		return 1;
	}

409
410
411
412
	list1 = asymmetric_indices(list1_raw, sym);
	list2 = asymmetric_indices(list2_raw, sym);

	/* Find common reflections and calculate ratio */
Thomas White's avatar
Thomas White committed
413
	ratio = reflist_new();
414
	ncom = 0;
Thomas White's avatar
Thomas White committed
415
416
	for ( refl1 = first_refl(list1, &iter);
	      refl1 != NULL;
Thomas White's avatar
Thomas White committed
417
418
	      refl1 = next_refl(refl1, iter) )
	{
Thomas White's avatar
Thomas White committed
419
		signed int h, k, l;
420
		double val1, val2;
Thomas White's avatar
Thomas White committed
421
422
		Reflection *refl2;
		Reflection *tr;
423

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

426
427
		refl2 = find_refl(list2, h, k, l);
		if ( refl2 == NULL ) continue;
428

429
		ncom++;
Thomas White's avatar
Thomas White committed
430
431
432
433
434
435
436

		val1 = get_intensity(refl1);
		val2 = get_intensity(refl2);

		/* Add divided version to 'output' list */
		tr = add_refl(ratio, h, k, l);
		set_int(tr, val1/val2);
437
		set_redundancy(tr, 1);
438
	}
Thomas White's avatar
Thomas White committed
439

440
441
442
443
444
	if ( ratiofile != NULL ) {
		write_reflist(ratiofile, ratio, cell);
	}
	reflist_free(ratio);

445
446
	gsl_set_error_handler_off();

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
	       num_reflections(list1), num_reflections(list2), ncom);

450
	R1 = stat_r1_ignore(list1, list2, &scale_r1fi, config_unity);
Thomas White's avatar
Thomas White committed
451
452
	STATUS("R1(F) = %5.4f %% (scale=%5.2e)"
	       " (ignoring negative intensities)\n",
453
	       R1*100.0, scale_r1fi);
Thomas White's avatar
Thomas White committed
454

455
	R1 = stat_r1_zero(list1, list2, &scale_r1, config_unity);
Thomas White's avatar
Thomas White committed
456
457
	STATUS("R1(F) = %5.4f %% (scale=%5.2e)"
	       " (zeroing negative intensities)\n",
458
	       R1*100.0, scale_r1);
Thomas White's avatar
Thomas White committed
459

460
	R2 = stat_r2(list1, list2, &scale_r2, config_unity);
461
	STATUS("R2(I) = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale_r2);
Thomas White's avatar
Thomas White committed
462

463
	R1i = stat_r1_i(list1, list2, &scale_r1i, config_unity);
464
	STATUS("R1(I) = %5.4f %% (scale=%5.2e)\n", R1i*100.0, scale_r1i);
Thomas White's avatar
Thomas White committed
465

466
	Rdiff = stat_rdiff_ignore(list1, list2, &scale_rdig,
467
	                          config_unity);
Thomas White's avatar
Thomas White committed
468
469
	STATUS("Rint(F) = %5.4f %% (scale=%5.2e)"
	       " (ignoring negative intensities)\n",
470
	       Rdiff*100.0, scale_rdig);
Thomas White's avatar
Thomas White committed
471

472
	Rdiff = stat_rdiff_zero(list1, list2, &scale, config_unity);
Thomas White's avatar
Thomas White committed
473
474
	STATUS("Rint(F) = %5.4f %% (scale=%5.2e)"
	       " (zeroing negative intensities)\n",
Thomas White's avatar
Thomas White committed
475
476
	       Rdiff*100.0, scale);

477
	Rdiff = stat_rdiff_intensity(list1, list2, &scale_rintint,
478
	                             config_unity);
Thomas White's avatar
Thomas White committed
479
	STATUS("Rint(I) = %5.4f %% (scale=%5.2e)\n",
480
	       Rdiff*100.0, scale_rintint);
Thomas White's avatar
Thomas White committed
481

482
	pearson = stat_pearson_i(list1, list2);
Thomas White's avatar
Thomas White committed
483
484
	STATUS("Pearson r(I) = %5.4f\n", pearson);

485
	pearson = stat_pearson_f_ignore(list1, list2);
Thomas White's avatar
Thomas White committed
486
487
488
	STATUS("Pearson r(F) = %5.4f (ignoring negative intensities)\n",
	       pearson);

489
	pearson = stat_pearson_f_zero(list1, list2);
Thomas White's avatar
Thomas White committed
490
491
	STATUS("Pearson r(F) = %5.4f (zeroing negative intensities)\n",
	       pearson);
Thomas White's avatar
Thomas White committed
492

493
494
495
496
497
498
499
500
501
502
	switch ( config_shells ) {
		case R_SHELL_R1I : scale_for_shells = scale_r1i;  break;
		case R_SHELL_R1F : scale_for_shells = scale_r1;  break;
		case R_SHELL_RSPLIT : scale_for_shells = scale_rintint;  break;
		default : scale_for_shells = 0.0;
	}

	if ( config_shells != R_SHELL_NONE ) {
		plot_shells(list1, list2, scale_for_shells,
		            cell, rmin_fix, rmax_fix, config_shells);
Thomas White's avatar
Thomas White committed
503
504
	}

505
506
	return 0;
}