compare_hkl.c 11.6 KB
Newer Older
1
2
3
4
5
/*
 * compare_hkl.c
 *
 * Compare reflection lists
 *
Thomas White's avatar
Thomas White committed
6
 * Copyright © 2012 Thomas White <taw@physics.org>
7
 *
Thomas White's avatar
Thomas White committed
8
9
10
11
12
13
14
15
16
17
18
19
20
21
 * 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/>.
22
23
24
25
26
27
28
29
30
31
32
33
34
35
 *
 */


#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
36
#include <assert.h>
37
#include <gsl/gsl_errno.h>
38
39

#include "utils.h"
Thomas White's avatar
Thomas White committed
40
#include "statistics.h"
41
#include "symmetry.h"
Thomas White's avatar
Thomas White committed
42
#include "reflist-utils.h"
43
44


Thomas White's avatar
Thomas White committed
45
/* Number of bins for plot of resolution shells */
46
#define NBINS (10)
Thomas White's avatar
Thomas White committed
47
48


49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
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);
}


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


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

	if ( cell == NULL ) {
Thomas White's avatar
Thomas White committed
109
		ERROR("Need the unit cell to plot resolution shells.\n");
Thomas White's avatar
Thomas White committed
110
111
112
113
114
		return;
	}

	for ( i=0; i<NBINS; i++ ) {
		num[i] = 0.0;
Thomas White's avatar
Thomas White committed
115
		cts[i] = 0;
116
117
		measured[i] = 0;
		measurements[i] = 0;
118
		snr[i] = 0;
Thomas White's avatar
Thomas White committed
119
120
	}

121
122
	/* Find resolution limits */
	resolution_limits(list1, cell, &rmin, &rmax);
123
	STATUS("1/d goes from %f to %f nm^-1\n", rmin/1e9, rmax/1e9);
124

125
126
	/* Widen the range just a little bit */
	rmin -= 0.001e9;
127
128
	rmax += 0.001e9;

129
130
131
	/* Fixed resolution shells if needed */
	if ( rmin_fix > 0.0 ) rmin = rmin_fix;
	if ( rmax_fix > 0.0 ) rmax = rmax_fix;
132

Thomas White's avatar
Thomas White committed
133
	/* Calculate the resolution bins */
134
135
136
137
138
139
140
141
142
143
	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);

144
		/* Shells of constant volume */
145
146
147
		rmaxs[i-1] = r;
		rmins[i] = r;

148
149
150
151
		/* Shells of constant thickness */
		//rmins[i] = rmins[i-1] + (rmax-rmin)/NBINS;
		//rmaxs[i-1] = rmins[i-1] + (rmax-rmin)/NBINS;

152
153
154
155
156
157
158
159
		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);

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

Thomas White's avatar
Thomas White committed
172
		get_indices(refl1, &h, &k, &l);
173
		d = 2.0 * resolution(cell, h, k, l);
Thomas White's avatar
Thomas White committed
174

175
176
177
178
179
180
181
		bin = -1;
		for ( j=0; j<NBINS; j++ ) {
			if ( (d>rmins[j]) && (d<=rmaxs[j]) ) {
				bin = j;
				break;
			}
		}
182
183

		/* Outside resolution range? */
184
185
186
187
		if ( bin == -1 ) {
			nout++;
			continue;
		}
Thomas White's avatar
Thomas White committed
188

189
190
191
		refl2 = find_refl(list2, h, k, l);
		if ( refl2 == NULL ) continue;

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

Thomas White's avatar
Thomas White committed
218
		ctot++;
219
		cts[bin]++;
Thomas White's avatar
Thomas White committed
220
221
	}

222
223
224
225
226
	if ( nout ) {
		STATUS("Warning; %i reflections outside resolution range.\n",
		       nout);
	}

227
228
229
230
231
232
	fh = fopen("shells.dat", "w");
	if ( fh == NULL ) {
		ERROR("Couldn't open 'shells.dat'\n");
		return;
	}

233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
	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;

	}
252

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

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

		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;

		}

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

	}

	fclose(fh);
Thomas White's avatar
Thomas White committed
281
282

	STATUS("Resolution shell information written to shells.dat.\n");
Thomas White's avatar
Thomas White committed
283
284
285
}


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

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

	/* Short options */
325
326
327
	while ((c = getopt_long(argc, argv, "ho:y:p:u",
	                        longopts, NULL)) != -1)
	{
328
329

		switch (c) {
Thomas White's avatar
Thomas White committed
330
		case 'h' :
331
332
333
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
334
		case 'o' :
Thomas White's avatar
Thomas White committed
335
			ratiofile = strdup(optarg);
336
337
			break;

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

342
343
344
345
		case 'p' :
			pdb = strdup(optarg);
			break;

346
347
348
349
		case 'u' :
			config_unity = 1;
			break;

Thomas White's avatar
Thomas White committed
350
		case 0 :
351
352
			break;

353
354
355
356
357
358
359
360
361
362
363
364
365
366
		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;

367
368
369
370
		case 4 :
			config_shells = get_r_shell(optarg);
			break;

Thomas White's avatar
Thomas White committed
371
		default :
372
373
374
375
376
			return 1;
		}

	}

377
378
379
380
381
	if ( argc != (optind+2) ) {
		ERROR("Please provide exactly two HKL files to compare.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
382
383
	if ( sym_str == NULL ) {
		sym_str = strdup("1");
384
	}
Thomas White's avatar
Thomas White committed
385
386
	sym = get_pointgroup(sym_str);
	free(sym_str);
387

388
389
390
	afile = strdup(argv[optind++]);
	bfile = strdup(argv[optind]);

391
392
393
394
395
396
397
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

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

398
399
	list1_raw = read_reflections(afile);
	if ( list1_raw == NULL ) {
Thomas White's avatar
Thomas White committed
400
		ERROR("Couldn't read file '%s'\n", afile);
401
402
		return 1;
	}
Thomas White's avatar
Thomas White committed
403

404
405
	list2_raw = read_reflections(bfile);
	if ( list2_raw == NULL ) {
Thomas White's avatar
Thomas White committed
406
		ERROR("Couldn't read file '%s'\n", bfile);
407
408
409
		return 1;
	}

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

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

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

439
440
		refl2 = find_refl(list2, h, k, l);
		if ( refl2 == NULL ) continue;
441

442
		ncom++;
Thomas White's avatar
Thomas White committed
443
444
445
446
447
448
449

		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);
450
		set_redundancy(tr, 1);
451
	}
Thomas White's avatar
Thomas White committed
452

453
454
455
456
457
	if ( ratiofile != NULL ) {
		write_reflist(ratiofile, ratio, cell);
	}
	reflist_free(ratio);

458
459
	gsl_set_error_handler_off();

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

463
	R1 = stat_r1_ignore(list1, list2, &scale_r1fi, config_unity);
Thomas White's avatar
Thomas White committed
464
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(list1, list2, &scale_r1, config_unity);
Thomas White's avatar
Thomas White committed
469
470
	STATUS("R1(F) = %5.4f %% (scale=%5.2e)"
	       " (zeroing negative intensities)\n",
471
	       R1*100.0, scale_r1);
Thomas White's avatar
Thomas White committed
472

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

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

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

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

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

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

498
	pearson = stat_pearson_f_ignore(list1, list2);
Thomas White's avatar
Thomas White committed
499
500
501
	STATUS("Pearson r(F) = %5.4f (ignoring negative intensities)\n",
	       pearson);

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

506
507
508
509
510
511
512
513
514
515
	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
516
517
	}

518
519
	return 0;
}