compare_hkl.c 10.7 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
220
221
	fh = fopen("shells.dat", "w");
	if ( fh == NULL ) {
		ERROR("Couldn't open 'shells.dat'\n");
		return;
	}

	fprintf(fh, "1/d centre   Rsplit / %%\n");

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

224
		double r, cen;
225
		cen = rmins[i] + (rmaxs[i] - rmins[i])/2.0;
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243

		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;

		}

244
245
		fprintf(fh, "%10.3f %10.2f %10i\n",
		        cen*1.0e-9, r*100.0, cts[i]);
Thomas White's avatar
Thomas White committed
246
247
248
249

	}

	fclose(fh);
Thomas White's avatar
Thomas White committed
250
251

	STATUS("Resolution shell information written to shells.dat.\n");
Thomas White's avatar
Thomas White committed
252
253
254
}


255
256
257
int main(int argc, char *argv[])
{
	int c;
258
	UnitCell *cell;
Thomas White's avatar
Thomas White committed
259
	char *ratiofile = NULL;
260
261
	char *afile = NULL;
	char *bfile = NULL;
Thomas White's avatar
Thomas White committed
262
263
	char *sym_str = NULL;
	SymOpList *sym;
264
	double scale, scale_r2, scale_rdig, R1, R2, R1i, Rdiff, pearson;
265
	double scale_rintint, scale_r1i, scale_r1, scale_r1fi;
Thomas White's avatar
Thomas White committed
266
267
268
	int ncom;
	RefList *list1;
	RefList *list2;
269
270
	RefList *list1_raw;
	RefList *list2_raw;
Thomas White's avatar
Thomas White committed
271
	RefList *ratio;
272
	enum r_shell config_shells = R_SHELL_NONE;
273
	char *pdb = NULL;
274
275
	float rmin_fix = -1.0;
	float rmax_fix = -1.0;
Thomas White's avatar
Thomas White committed
276
277
	Reflection *refl1;
	RefListIterator *iter;
278
	int config_unity = 0;
279
	double scale_for_shells;
280
281
282
283

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
Thomas White's avatar
Thomas White committed
284
		{"ratio" ,             1, NULL,               'o'},
285
		{"symmetry",           1, NULL,               'y'},
286
		{"shells",             1, NULL,               4},
287
		{"pdb",                1, NULL,               'p'},
288
289
		{"rmin",               1, NULL,               2},
		{"rmax",               1, NULL,               3},
290
291
292
293
		{0, 0, NULL, 0}
	};

	/* Short options */
294
295
296
	while ((c = getopt_long(argc, argv, "ho:y:p:u",
	                        longopts, NULL)) != -1)
	{
297
298

		switch (c) {
Thomas White's avatar
Thomas White committed
299
		case 'h' :
300
301
302
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
303
		case 'o' :
Thomas White's avatar
Thomas White committed
304
			ratiofile = strdup(optarg);
305
306
			break;

307
		case 'y' :
Thomas White's avatar
Thomas White committed
308
			sym_str = strdup(optarg);
309
310
			break;

311
312
313
314
		case 'p' :
			pdb = strdup(optarg);
			break;

315
316
317
318
		case 'u' :
			config_unity = 1;
			break;

Thomas White's avatar
Thomas White committed
319
		case 0 :
320
321
			break;

322
323
324
325
326
327
328
329
330
331
332
333
334
335
		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;

336
337
338
339
		case 4 :
			config_shells = get_r_shell(optarg);
			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;
	}

Thomas White's avatar
Thomas White committed
351
352
	if ( sym_str == NULL ) {
		sym_str = strdup("1");
353
	}
Thomas White's avatar
Thomas White committed
354
355
	sym = get_pointgroup(sym_str);
	free(sym_str);
356

357
358
359
	afile = strdup(argv[optind++]);
	bfile = strdup(argv[optind]);

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

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

367
368
	list1_raw = read_reflections(afile);
	if ( list1_raw == NULL ) {
Thomas White's avatar
Thomas White committed
369
		ERROR("Couldn't read file '%s'\n", afile);
370
371
		return 1;
	}
Thomas White's avatar
Thomas White committed
372

373
374
	list2_raw = read_reflections(bfile);
	if ( list2_raw == NULL ) {
Thomas White's avatar
Thomas White committed
375
		ERROR("Couldn't read file '%s'\n", bfile);
376
377
378
		return 1;
	}

Thomas White's avatar
Thomas White committed
379
	/* Check that the intensities have the correct symmetry */
380
	if ( check_list_symmetry(list1_raw, sym) ) {
Thomas White's avatar
Thomas White committed
381
		ERROR("The first input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
382
		      " have symmetry %s\n", symmetry_name(sym));
Thomas White's avatar
Thomas White committed
383
384
		return 1;
	}
385
	if ( check_list_symmetry(list2_raw, sym) ) {
Thomas White's avatar
Thomas White committed
386
		ERROR("The second input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
387
		      " have symmetry %s\n", symmetry_name(sym));
Thomas White's avatar
Thomas White committed
388
389
390
		return 1;
	}

391
392
393
394
	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
395
	ratio = reflist_new();
396
	ncom = 0;
Thomas White's avatar
Thomas White committed
397
398
	for ( refl1 = first_refl(list1, &iter);
	      refl1 != NULL;
Thomas White's avatar
Thomas White committed
399
400
	      refl1 = next_refl(refl1, iter) )
	{
Thomas White's avatar
Thomas White committed
401
		signed int h, k, l;
402
		double val1, val2;
Thomas White's avatar
Thomas White committed
403
404
		Reflection *refl2;
		Reflection *tr;
405

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

408
409
		refl2 = find_refl(list2, h, k, l);
		if ( refl2 == NULL ) continue;
410

411
		ncom++;
Thomas White's avatar
Thomas White committed
412
413
414
415
416
417
418

		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);
419
		set_redundancy(tr, 1);
420
	}
Thomas White's avatar
Thomas White committed
421

422
423
424
425
426
	if ( ratiofile != NULL ) {
		write_reflist(ratiofile, ratio, cell);
	}
	reflist_free(ratio);

427
428
	gsl_set_error_handler_off();

Thomas White's avatar
Thomas White committed
429
	STATUS("%i,%i reflections: %i in common\n",
Thomas White's avatar
Thomas White committed
430
431
	       num_reflections(list1), num_reflections(list2), ncom);

432
	R1 = stat_r1_ignore(list1, list2, &scale_r1fi, config_unity);
Thomas White's avatar
Thomas White committed
433
434
	STATUS("R1(F) = %5.4f %% (scale=%5.2e)"
	       " (ignoring negative intensities)\n",
435
	       R1*100.0, scale_r1fi);
Thomas White's avatar
Thomas White committed
436

437
	R1 = stat_r1_zero(list1, list2, &scale_r1, config_unity);
Thomas White's avatar
Thomas White committed
438
439
	STATUS("R1(F) = %5.4f %% (scale=%5.2e)"
	       " (zeroing negative intensities)\n",
440
	       R1*100.0, scale_r1);
Thomas White's avatar
Thomas White committed
441

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

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

448
	Rdiff = stat_rdiff_ignore(list1, list2, &scale_rdig,
449
	                          config_unity);
Thomas White's avatar
Thomas White committed
450
451
	STATUS("Rint(F) = %5.4f %% (scale=%5.2e)"
	       " (ignoring negative intensities)\n",
452
	       Rdiff*100.0, scale_rdig);
Thomas White's avatar
Thomas White committed
453

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

459
	Rdiff = stat_rdiff_intensity(list1, list2, &scale_rintint,
460
	                             config_unity);
Thomas White's avatar
Thomas White committed
461
	STATUS("Rint(I) = %5.4f %% (scale=%5.2e)\n",
462
	       Rdiff*100.0, scale_rintint);
Thomas White's avatar
Thomas White committed
463

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

467
	pearson = stat_pearson_f_ignore(list1, list2);
Thomas White's avatar
Thomas White committed
468
469
470
	STATUS("Pearson r(F) = %5.4f (ignoring negative intensities)\n",
	       pearson);

471
	pearson = stat_pearson_f_zero(list1, list2);
Thomas White's avatar
Thomas White committed
472
473
	STATUS("Pearson r(F) = %5.4f (zeroing negative intensities)\n",
	       pearson);
Thomas White's avatar
Thomas White committed
474

475
476
477
478
479
480
481
482
483
484
	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
485
486
	}

487
488
	return 0;
}