compare_hkl.c 7.02 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 */
Thomas White's avatar
Thomas White committed
32
33
34
#define NBINS (10)


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
51
static void plot_shells(const double *ref1, const double *ref2,
                        ReflItemList *items, double scale, UnitCell *cell)
Thomas White's avatar
Thomas White committed
52
53
54
55
56
57
58
59
{
	double num[NBINS];
	double den[NBINS];
	double rmin, rmax, rstep;
	int i;
	FILE *fh;

	if ( cell == NULL ) {
Thomas White's avatar
Thomas White committed
60
		ERROR("Need the unit cell to plot resolution shells.\n");
Thomas White's avatar
Thomas White committed
61
62
63
		return;
	}

Thomas White's avatar
Thomas White committed
64
	fh = fopen("shells.dat", "w");
Thomas White's avatar
Thomas White committed
65
	if ( fh == NULL ) {
Thomas White's avatar
Thomas White committed
66
		ERROR("Couldn't open 'shells.dat'\n");
Thomas White's avatar
Thomas White committed
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
		return;
	}

	for ( i=0; i<NBINS; i++ ) {
		num[i] = 0.0;
		den[i] = 0.0;
	}

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

		struct refl_item *it;
		signed int h, k, l;
		double res;

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

		res = 2.0*resolution(cell, h, k, l);
		if ( res > rmax ) rmax = res;
		if ( res < rmin ) rmin = res;

	}
	rstep = (rmax-rmin) / NBINS;

	for ( i=0; i<num_items(items); i++ ) {

		struct refl_item *it;
		signed int h, k, l;
		double res;
		int bin;
		double i1, i2;

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

		res = 2.0*resolution(cell, h, k, l);

		bin = (res-rmin)/rstep;

		i1 = lookup_intensity(ref1, h, k, l);
		i2 = scale * lookup_intensity(ref2, h, k, l);

		num[bin] += pow(i1 - i2, 2.0);
		den[bin] += pow(i1, 2.0);

	}

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

		double r2, cen;
		cen = rmin + rstep*i + rstep/2.0;
		r2 = sqrt(num[i]/den[i]);
		fprintf(fh, "%f %f\n", cen/1.0e9, r2*100.0);

	}

	fclose(fh);
}


129
130
131
132
133
int main(int argc, char *argv[])
{
	int c;
	double *ref1;
	double *ref2;
134
	double *ref2_transformed;
135
	double *out;
136
	UnitCell *cell;
137
138
139
	char *outfile = NULL;
	char *afile = NULL;
	char *bfile = NULL;
Thomas White's avatar
Thomas White committed
140
	char *sym = NULL;
141
	double scale, scale_r2, R1, R2, R1i, Rdiff, pearson;
Thomas White's avatar
Thomas White committed
142
143
	int i, ncom;
	ReflItemList *i1, *i2, *icommon;
Thomas White's avatar
Thomas White committed
144
	int config_shells = 0;
145
	char *pdb = NULL;
146
147
148
149
	double *esd1;
	double *esd2;
	int rej1 = 0;
	int rej2 = 0;
150
151
152
153
154

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"output",             1, NULL,               'o'},
155
		{"symmetry",           1, NULL,               'y'},
Thomas White's avatar
Thomas White committed
156
		{"shells",             0, &config_shells,     1},
157
		{"pdb",                1, NULL,               'p'},
158
159
160
161
		{0, 0, NULL, 0}
	};

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

		switch (c) {
Thomas White's avatar
Thomas White committed
165
		case 'h' :
166
167
168
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
169
		case 'o' :
170
171
172
			outfile = strdup(optarg);
			break;

173
174
175
176
		case 'y' :
			sym = strdup(optarg);
			break;

177
178
179
180
		case 'p' :
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
181
		case 0 :
182
183
			break;

Thomas White's avatar
Thomas White committed
184
		default :
185
186
187
188
189
			return 1;
		}

	}

190
191
192
193
194
	if ( argc != (optind+2) ) {
		ERROR("Please provide exactly two HKL files to compare.\n");
		return 1;
	}

195
196
197
198
	if ( sym == NULL ) {
		sym = strdup("1");
	}

199
200
201
	afile = strdup(argv[optind++]);
	bfile = strdup(argv[optind]);

202
203
204
205
206
207
208
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

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

Thomas White's avatar
Thomas White committed
209
	ref1 = new_list_intensity();
210
211
	esd1 = new_list_sigma();
	i1 = read_reflections(afile, ref1, NULL, NULL, esd1);
Thomas White's avatar
Thomas White committed
212
	if ( i1 == NULL ) {
213
214
215
		ERROR("Couldn't open file '%s'\n", afile);
		return 1;
	}
Thomas White's avatar
Thomas White committed
216
	ref2 = new_list_intensity();
217
218
	esd2 = new_list_sigma();
	i2 = read_reflections(bfile, ref2, NULL, NULL, esd2);
Thomas White's avatar
Thomas White committed
219
	if ( i2 == NULL ) {
220
221
222
223
		ERROR("Couldn't open file '%s'\n", bfile);
		return 1;
	}

Thomas White's avatar
Thomas White committed
224
225
	/* List for output scale factor map */
	out = new_list_intensity();
226

227
228
229
230
231
	/* 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
232
233
		struct refl_item *it;
		signed int h, k, l;
234
235
		signed int he, ke, le;
		double val1, val2;
236
237
		double sig1, sig2;
		int ig = 0;
238

239
		it = get_item(i1, i);
Thomas White's avatar
Thomas White committed
240
		h = it->h;  k = it->k;  l = it->l;
241

242
		if ( !find_unique_equiv(i2, h, k, l, sym, &he, &ke, &le) ) {
Thomas White's avatar
Thomas White committed
243
244
			//STATUS("%i %i %i not matched (%f nm).\n", h, k, l,
			//       1.0/(2.0*resolution(cell, h, k, l)/1e9));
245
246
			continue;
		}
247

248
249
		val1 = lookup_intensity(ref1, h, k, l);
		val2 = lookup_intensity(ref2, he, ke, le);
250
		sig1 = lookup_sigma(esd1, h, k, l);
251
		sig2 = lookup_sigma(esd2, he, ke, le);
252
253
254
255
256
257
258
259
260
261
262

		if ( val1 < 3.0 * sig1 ) {
			rej1++;
			ig = 1;
		}
		if ( val2 < 3.0 * sig2 ) {
			rej2++;
			ig = 1;
		}
		if ( ig ) continue;

263
264
265
		set_intensity(ref2_transformed, h, k, l, val2);
		set_intensity(out, h, k, l, val1/val2);
		add_item(icommon, h, k, l);
266
267

	}
268
	ncom = num_items(icommon);
269
270
271
272
	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);
273

Thomas White's avatar
Thomas White committed
274
275
	STATUS("%i,%i reflections: %i in common\n",
	       num_items(i1), num_items(i2), ncom);
Thomas White's avatar
Thomas White committed
276
277

	R1 = stat_r1_ignore(ref1, ref2_transformed, icommon, &scale);
278
	STATUS("R1(F) = %5.4f %% (scale=%5.2e) (ignoring negative intensities)\n",
Thomas White's avatar
Thomas White committed
279
280
281
	       R1*100.0, scale);

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

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

288
289
	R1i = stat_r1_i(ref1, ref2_transformed, icommon, &scale);
	STATUS("R1(I) = %5.4f %% (scale=%5.2e)\n", R1i*100.0, scale);
Thomas White's avatar
Thomas White committed
290
291

	Rdiff = stat_rdiff_ignore(ref1, ref2_transformed, icommon, &scale);
292
	STATUS("Rdiff(F) = %5.4f %% (scale=%5.2e) (ignoring negative intensities)\n",
Thomas White's avatar
Thomas White committed
293
294
295
	       Rdiff*100.0, scale);

	Rdiff = stat_rdiff_zero(ref1, ref2_transformed, icommon, &scale);
296
	STATUS("Rdiff(F) = %5.4f %% (scale=%5.2e) (zeroing negative intensities)\n",
Thomas White's avatar
Thomas White committed
297
298
299
300
301
302
303
304
305
306
307
308
	       Rdiff*100.0, scale);

	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
309

Thomas White's avatar
Thomas White committed
310
311
	if ( config_shells ) {
		plot_shells(ref1, ref2_transformed, icommon, scale_r2, cell);
Thomas White's avatar
Thomas White committed
312
313
	}

314
	if ( outfile != NULL ) {
Thomas White's avatar
Thomas White committed
315
		write_reflections(outfile, icommon, out, NULL, NULL, cell);
316
	}
317
318
319

	return 0;
}