compare_hkl.c 11.7 KB
Newer Older
1
2
3
4
5
/*
 * compare_hkl.c
 *
 * Compare reflection lists
 *
6
7
8
9
10
 * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
 *                  a research centre of the Helmholtz Association.
 *
 * Authors:
 *   2009-2012 Thomas White <taw@physics.org>
11
 *
Thomas White's avatar
Thomas White committed
12
13
14
15
16
17
18
19
20
21
22
23
24
25
 * This file is part of CrystFEL.
 *
 * CrystFEL is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * CrystFEL is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
26
27
28
29
30
31
32
33
34
35
36
37
38
39
 *
 */


#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
40
#include <assert.h>
41
#include <gsl/gsl_errno.h>
42
43

#include "utils.h"
Thomas White's avatar
Thomas White committed
44
#include "statistics.h"
45
#include "symmetry.h"
Thomas White's avatar
Thomas White committed
46
#include "reflist-utils.h"
47
48


Thomas White's avatar
Thomas White committed
49
/* Number of bins for plot of resolution shells */
50
#define NBINS (10)
Thomas White's avatar
Thomas White committed
51
52


53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
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);
}


74
75
static void show_help(const char *s)
{
76
	printf("Syntax: %s [options] <file1.hkl> <file2.hkl>\n\n", s);
77
78
79
80
	printf(
"Compare intensity lists.\n"
"\n"
"  -h, --help                 Display this help message.\n"
Thomas White's avatar
Thomas White committed
81
"  -o, --ratio=<filename>     Specify output filename for ratios.\n"
82
"  -y, --symmetry=<sym>       The symmetry of both the input files.\n"
Thomas White's avatar
Thomas White committed
83
"  -p, --pdb=<filename>       PDB file to use.\n"
84
85
"      --shells=<FoM>         Plot this figure of merit in resolution shells.\n"
"                              Choose from: 'Rsplit', 'R1f' and 'R1i'.\n"
86
87
"      --rmin=<res>           Fix lower resolution limit for --shells (m^-1).\n"
"      --rmax=<res>           Fix upper resolution limit for --shells (m^-1).\n"
88
89
90
91
"\n");
}


92
static void plot_shells(RefList *list1, RefList *list2, double scale,
93
94
                        UnitCell *cell, double rmin_fix, double rmax_fix,
                        enum r_shell config_shells)
Thomas White's avatar
Thomas White committed
95
96
{
	double num[NBINS];
Thomas White's avatar
Thomas White committed
97
	int cts[NBINS];
98
99
100
101
102
	unsigned int measurements[NBINS];
	unsigned int measured[NBINS];
	double total_vol, vol_per_shell;
	double rmins[NBINS];
	double rmaxs[NBINS];
103
	double snr[NBINS];
104
	double rmin, rmax;
Thomas White's avatar
Thomas White committed
105
	int i;
Thomas White's avatar
Thomas White committed
106
107
	Reflection *refl1;
	RefListIterator *iter;
108
109
110
	FILE *fh;
	double den;
	int ctot, nout;
Thomas White's avatar
Thomas White committed
111
112

	if ( cell == NULL ) {
Thomas White's avatar
Thomas White committed
113
		ERROR("Need the unit cell to plot resolution shells.\n");
Thomas White's avatar
Thomas White committed
114
115
116
117
118
		return;
	}

	for ( i=0; i<NBINS; i++ ) {
		num[i] = 0.0;
Thomas White's avatar
Thomas White committed
119
		cts[i] = 0;
120
121
		measured[i] = 0;
		measurements[i] = 0;
122
		snr[i] = 0;
Thomas White's avatar
Thomas White committed
123
124
	}

125
126
	/* Find resolution limits */
	resolution_limits(list1, cell, &rmin, &rmax);
127
	STATUS("1/d goes from %f to %f nm^-1\n", rmin/1e9, rmax/1e9);
128

129
130
	/* Widen the range just a little bit */
	rmin -= 0.001e9;
131
132
	rmax += 0.001e9;

133
134
135
	/* Fixed resolution shells if needed */
	if ( rmin_fix > 0.0 ) rmin = rmin_fix;
	if ( rmax_fix > 0.0 ) rmax = rmax_fix;
136

Thomas White's avatar
Thomas White committed
137
	/* Calculate the resolution bins */
138
139
140
141
142
143
144
145
146
147
	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);

148
		/* Shells of constant volume */
149
150
151
		rmaxs[i-1] = r;
		rmins[i] = r;

152
153
154
155
		/* Shells of constant thickness */
		//rmins[i] = rmins[i-1] + (rmax-rmin)/NBINS;
		//rmaxs[i-1] = rmins[i-1] + (rmax-rmin)/NBINS;

156
157
158
159
160
161
162
163
		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);

164
	den = 0.0;  ctot = 0;  nout = 0;
Thomas White's avatar
Thomas White committed
165
166
	for ( refl1 = first_refl(list1, &iter);
	      refl1 != NULL;
167
168
	      refl1 = next_refl(refl1, iter) )
	{
Thomas White's avatar
Thomas White committed
169
		signed int h, k, l;
170
		double d;
Thomas White's avatar
Thomas White committed
171
		int bin;
172
		double i1, i2, f1, f2;
173
		int j;
174
		Reflection *refl2;
Thomas White's avatar
Thomas White committed
175

Thomas White's avatar
Thomas White committed
176
		get_indices(refl1, &h, &k, &l);
177
		d = 2.0 * resolution(cell, h, k, l);
Thomas White's avatar
Thomas White committed
178

179
180
181
182
183
184
185
		bin = -1;
		for ( j=0; j<NBINS; j++ ) {
			if ( (d>rmins[j]) && (d<=rmaxs[j]) ) {
				bin = j;
				break;
			}
		}
186
187

		/* Outside resolution range? */
188
189
190
191
		if ( bin == -1 ) {
			nout++;
			continue;
		}
Thomas White's avatar
Thomas White committed
192

193
194
195
		refl2 = find_refl(list2, h, k, l);
		if ( refl2 == NULL ) continue;

Thomas White's avatar
Thomas White committed
196
		i1 = get_intensity(refl1);
197
198
199
200
201
202
203
204
		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);
205
			den += i1 + scale*i2;
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
			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
221

Thomas White's avatar
Thomas White committed
222
		ctot++;
223
		cts[bin]++;
Thomas White's avatar
Thomas White committed
224
225
	}

226
227
228
229
230
	if ( nout ) {
		STATUS("Warning; %i reflections outside resolution range.\n",
		       nout);
	}

231
232
233
234
235
236
	fh = fopen("shells.dat", "w");
	if ( fh == NULL ) {
		ERROR("Couldn't open 'shells.dat'\n");
		return;
	}

237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
	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;

	}
256

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

259
		double r, cen;
260
		cen = rmins[i] + (rmaxs[i] - rmins[i])/2.0;
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278

		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;

		}

279
280
		fprintf(fh, "%10.3f %10.2f %10i\n",
		        cen*1.0e-9, r*100.0, cts[i]);
Thomas White's avatar
Thomas White committed
281
282
283
284

	}

	fclose(fh);
Thomas White's avatar
Thomas White committed
285
286

	STATUS("Resolution shell information written to shells.dat.\n");
Thomas White's avatar
Thomas White committed
287
288
289
}


290
291
292
int main(int argc, char *argv[])
{
	int c;
293
	UnitCell *cell;
Thomas White's avatar
Thomas White committed
294
	char *ratiofile = NULL;
295
296
	char *afile = NULL;
	char *bfile = NULL;
Thomas White's avatar
Thomas White committed
297
298
	char *sym_str = NULL;
	SymOpList *sym;
299
	double scale, scale_r2, scale_rdig, R1, R2, R1i, Rdiff, pearson;
300
	double scale_rintint, scale_r1i, scale_r1, scale_r1fi;
Thomas White's avatar
Thomas White committed
301
302
303
	int ncom;
	RefList *list1;
	RefList *list2;
304
305
	RefList *list1_raw;
	RefList *list2_raw;
Thomas White's avatar
Thomas White committed
306
	RefList *ratio;
307
	enum r_shell config_shells = R_SHELL_NONE;
308
	char *pdb = NULL;
309
310
	float rmin_fix = -1.0;
	float rmax_fix = -1.0;
Thomas White's avatar
Thomas White committed
311
312
	Reflection *refl1;
	RefListIterator *iter;
313
	int config_unity = 0;
314
	double scale_for_shells;
315
316
317
318

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
Thomas White's avatar
Thomas White committed
319
		{"ratio" ,             1, NULL,               'o'},
320
		{"symmetry",           1, NULL,               'y'},
321
		{"shells",             1, NULL,               4},
322
		{"pdb",                1, NULL,               'p'},
323
324
		{"rmin",               1, NULL,               2},
		{"rmax",               1, NULL,               3},
325
326
327
328
		{0, 0, NULL, 0}
	};

	/* Short options */
329
330
331
	while ((c = getopt_long(argc, argv, "ho:y:p:u",
	                        longopts, NULL)) != -1)
	{
332
333

		switch (c) {
Thomas White's avatar
Thomas White committed
334
		case 'h' :
335
336
337
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
338
		case 'o' :
Thomas White's avatar
Thomas White committed
339
			ratiofile = strdup(optarg);
340
341
			break;

342
		case 'y' :
Thomas White's avatar
Thomas White committed
343
			sym_str = strdup(optarg);
344
345
			break;

346
347
348
349
		case 'p' :
			pdb = strdup(optarg);
			break;

350
351
352
353
		case 'u' :
			config_unity = 1;
			break;

Thomas White's avatar
Thomas White committed
354
		case 0 :
355
356
			break;

357
358
359
360
361
362
363
364
365
366
367
368
369
370
		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;

371
372
373
374
		case 4 :
			config_shells = get_r_shell(optarg);
			break;

Thomas White's avatar
Thomas White committed
375
		default :
376
377
378
379
380
			return 1;
		}

	}

381
382
383
384
385
	if ( argc != (optind+2) ) {
		ERROR("Please provide exactly two HKL files to compare.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
386
387
	if ( sym_str == NULL ) {
		sym_str = strdup("1");
388
	}
Thomas White's avatar
Thomas White committed
389
390
	sym = get_pointgroup(sym_str);
	free(sym_str);
391

392
393
394
	afile = strdup(argv[optind++]);
	bfile = strdup(argv[optind]);

395
396
397
398
399
400
401
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

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

402
403
	list1_raw = read_reflections(afile);
	if ( list1_raw == NULL ) {
Thomas White's avatar
Thomas White committed
404
		ERROR("Couldn't read file '%s'\n", afile);
405
406
		return 1;
	}
Thomas White's avatar
Thomas White committed
407

408
409
	list2_raw = read_reflections(bfile);
	if ( list2_raw == NULL ) {
Thomas White's avatar
Thomas White committed
410
		ERROR("Couldn't read file '%s'\n", bfile);
411
412
413
		return 1;
	}

Thomas White's avatar
Thomas White committed
414
	/* Check that the intensities have the correct symmetry */
415
	if ( check_list_symmetry(list1_raw, sym) ) {
Thomas White's avatar
Thomas White committed
416
		ERROR("The first input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
417
		      " have symmetry %s\n", symmetry_name(sym));
Thomas White's avatar
Thomas White committed
418
419
		return 1;
	}
420
	if ( check_list_symmetry(list2_raw, sym) ) {
Thomas White's avatar
Thomas White committed
421
		ERROR("The second input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
422
		      " have symmetry %s\n", symmetry_name(sym));
Thomas White's avatar
Thomas White committed
423
424
425
		return 1;
	}

426
427
428
429
	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
430
	ratio = reflist_new();
431
	ncom = 0;
Thomas White's avatar
Thomas White committed
432
433
	for ( refl1 = first_refl(list1, &iter);
	      refl1 != NULL;
Thomas White's avatar
Thomas White committed
434
435
	      refl1 = next_refl(refl1, iter) )
	{
Thomas White's avatar
Thomas White committed
436
		signed int h, k, l;
437
		double val1, val2;
Thomas White's avatar
Thomas White committed
438
439
		Reflection *refl2;
		Reflection *tr;
440

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

443
444
		refl2 = find_refl(list2, h, k, l);
		if ( refl2 == NULL ) continue;
445

446
		ncom++;
Thomas White's avatar
Thomas White committed
447
448
449
450
451
452

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

		/* Add divided version to 'output' list */
		tr = add_refl(ratio, h, k, l);
453
		set_intensity(tr, val1/val2);
454
		set_redundancy(tr, 1);
455
	}
Thomas White's avatar
Thomas White committed
456

457
	if ( ratiofile != NULL ) {
458
		write_reflist(ratiofile, ratio);
459
460
461
	}
	reflist_free(ratio);

462
463
	gsl_set_error_handler_off();

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

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

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

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

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

483
	Rdiff = stat_rdiff_ignore(list1, list2, &scale_rdig,
484
	                          config_unity);
Thomas White's avatar
Thomas White committed
485
486
	STATUS("Rint(F) = %5.4f %% (scale=%5.2e)"
	       " (ignoring negative intensities)\n",
487
	       Rdiff*100.0, scale_rdig);
Thomas White's avatar
Thomas White committed
488

489
	Rdiff = stat_rdiff_zero(list1, list2, &scale, config_unity);
Thomas White's avatar
Thomas White committed
490
491
	STATUS("Rint(F) = %5.4f %% (scale=%5.2e)"
	       " (zeroing negative intensities)\n",
Thomas White's avatar
Thomas White committed
492
493
	       Rdiff*100.0, scale);

494
	Rdiff = stat_rdiff_intensity(list1, list2, &scale_rintint,
495
	                             config_unity);
Thomas White's avatar
Thomas White committed
496
	STATUS("Rint(I) = %5.4f %% (scale=%5.2e)\n",
497
	       Rdiff*100.0, scale_rintint);
Thomas White's avatar
Thomas White committed
498

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

502
	pearson = stat_pearson_f_ignore(list1, list2);
Thomas White's avatar
Thomas White committed
503
504
505
	STATUS("Pearson r(F) = %5.4f (ignoring negative intensities)\n",
	       pearson);

506
	pearson = stat_pearson_f_zero(list1, list2);
Thomas White's avatar
Thomas White committed
507
508
	STATUS("Pearson r(F) = %5.4f (zeroing negative intensities)\n",
	       pearson);
Thomas White's avatar
Thomas White committed
509

510
511
512
513
514
515
516
517
518
519
	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
520
521
	}

522
523
	return 0;
}