compare_hkl.c 5.18 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
32
33
34
/* Number of bins for Luzzati plot */
#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
45
46
47
"\n");
}


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

	if ( cell == NULL ) {
		ERROR("Need the unit cell to plot the Luzzati plot.\n");
		return;
	}

	fh = fopen("luzzati.dat", "w");
	if ( fh == NULL ) {
		ERROR("Couldn't open 'luzzati.dat'\n");
		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);
}


127
128
129
130
131
int main(int argc, char *argv[])
{
	int c;
	double *ref1;
	double *ref2;
132
	double *ref2_transformed;
133
	double *out;
134
	UnitCell *cell;
135
136
137
	char *outfile = NULL;
	char *afile = NULL;
	char *bfile = NULL;
Thomas White's avatar
Thomas White committed
138
	char *sym = NULL;
Thomas White's avatar
Thomas White committed
139
	double scale, scale_r2, R1, R2, Rdiff, pearson;
Thomas White's avatar
Thomas White committed
140
141
	int i, ncom;
	ReflItemList *i1, *i2, *icommon;
Thomas White's avatar
Thomas White committed
142
	int config_luzzati = 0;
143
144
145
146
147

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"output",             1, NULL,               'o'},
148
		{"symmetry",           1, NULL,               'y'},
Thomas White's avatar
Thomas White committed
149
		{"luzzati",            0, &config_luzzati,     1},
150
151
152
153
		{0, 0, NULL, 0}
	};

	/* Short options */
154
	while ((c = getopt_long(argc, argv, "ho:y:", longopts, NULL)) != -1) {
155
156

		switch (c) {
Thomas White's avatar
Thomas White committed
157
		case 'h' :
158
159
160
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
161
		case 'o' :
162
163
164
			outfile = strdup(optarg);
			break;

165
166
167
168
		case 'y' :
			sym = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
169
		case 0 :
170
171
			break;

Thomas White's avatar
Thomas White committed
172
		default :
173
174
175
176
177
			return 1;
		}

	}

178
179
180
181
182
	if ( argc != (optind+2) ) {
		ERROR("Please provide exactly two HKL files to compare.\n");
		return 1;
	}

183
184
185
186
	if ( sym == NULL ) {
		sym = strdup("1");
	}

187
188
189
	afile = strdup(argv[optind++]);
	bfile = strdup(argv[optind]);

190
	cell = load_cell_from_pdb("molecule.pdb");
Thomas White's avatar
Thomas White committed
191
192
	ref1 = new_list_intensity();
	i1 = read_reflections(afile, ref1, NULL, NULL);
Thomas White's avatar
Thomas White committed
193
	if ( i1 == NULL ) {
194
195
196
		ERROR("Couldn't open file '%s'\n", afile);
		return 1;
	}
Thomas White's avatar
Thomas White committed
197
198
	ref2 = new_list_intensity();
	i2 = read_reflections(bfile, ref2, NULL, NULL);
Thomas White's avatar
Thomas White committed
199
	if ( i2 == NULL ) {
200
201
202
203
		ERROR("Couldn't open file '%s'\n", bfile);
		return 1;
	}

Thomas White's avatar
Thomas White committed
204
205
	/* List for output scale factor map */
	out = new_list_intensity();
206

207
208
209
210
211
	/* 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
212
213
		struct refl_item *it;
		signed int h, k, l;
214
215
		signed int he, ke, le;
		double val1, val2;
216

217
		it = get_item(i1, i);
Thomas White's avatar
Thomas White committed
218
		h = it->h;  k = it->k;  l = it->l;
219

220
221
222
		if ( !find_unique_equiv(i2, h, k, l, sym, &he, &ke, &le) ) {
			continue;
		}
223

224
225
226
227
228
		val1 = lookup_intensity(ref1, h, k, l);
		val2 = lookup_intensity(ref2, he, ke, le);
		set_intensity(ref2_transformed, h, k, l, val2);
		set_intensity(out, h, k, l, val1/val2);
		add_item(icommon, h, k, l);
229
230

	}
231
	ncom = num_items(icommon);
232

Thomas White's avatar
Thomas White committed
233
234
	STATUS("%i,%i reflections: %i in common\n",
	       num_items(i1), num_items(i2), ncom);
Thomas White's avatar
Thomas White committed
235
236
	R1 = stat_r1(ref1, ref2_transformed, icommon, &scale);
	STATUS("R1 = %5.4f %% (scale=%5.2e)\n", R1*100.0, scale);
Thomas White's avatar
Thomas White committed
237
238
	R2 = stat_r2(ref1, ref2_transformed, icommon, &scale_r2);
	STATUS("R2 = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale_r2);
Thomas White's avatar
Thomas White committed
239
240
	Rdiff = stat_rdiff(ref1, ref2_transformed, icommon, &scale);
	STATUS("Rdiff = %5.4f %% (scale=%5.2e)\n", Rdiff*100.0, scale);
241
	pearson = stat_pearson(ref1, ref2_transformed, icommon);
242
	STATUS("Pearson r = %5.4f\n", pearson);
Thomas White's avatar
Thomas White committed
243

Thomas White's avatar
Thomas White committed
244
245
246
247
	if ( config_luzzati ) {
		plot_luzzati(ref1, ref2_transformed, icommon, scale_r2, cell);
	}

248
	if ( outfile != NULL ) {
Thomas White's avatar
Thomas White committed
249
		write_reflections(outfile, icommon, out, NULL, NULL, cell);
250
	}
251
252
253

	return 0;
}