process_hkl.c 10.4 KB
Newer Older
1
2
3
4
5
/*
 * process_hkl.c
 *
 * Assemble and process FEL Bragg intensities
 *
Thomas White's avatar
Thomas White committed
6
 * (c) 2006-2010 Thomas White <taw@physics.org>
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
 *
 * 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 "statistics.h"
#include "sfac.h"
27
#include "reflections.h"
28
#include "likelihood.h"
29
30


Thomas White's avatar
Thomas White committed
31
32
33
34
/* Number of divisions for R vs |q| graphs */
#define RVQDV (20)


35
36
37
static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
Thomas White's avatar
Thomas White committed
38
39
40
	printf(
"Assemble and process FEL Bragg intensities.\n"
"\n"
41
42
43
44
"  -h, --help                Display this help message.\n"
"  -i, --input=<filename>    Specify input filename (\"-\" for stdin).\n"
"  -o, --output=<filename>   Specify output filename for merged intensities\n"
"                             (don't specify for no output).\n"
Thomas White's avatar
Thomas White committed
45
"\n"
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
"      --max-only            Take the integrated intensity to be equal to the\n"
"                             maximum intensity measured for that reflection.\n"
"                             The default is to use the mean value from all\n"
"                             measurements.\n"
"      --sum                 Sum (rather than average) the intensities for the\n"
"                             final output list.  This is useful for comparing\n"
"                             results to radially summed powder patterns, but\n"
"                             will break R-factor analysis.\n"
"      --stop-after=<n>      Stop after processing n patterns.  Zero means\n"
"                             keep going until the end of the input, and is\n"
"                             the default.\n"
"  -c, --compare-with=<file> Compare with reflection intensities in this file\n"
"\n"
"  -e, --output-every=<n>    Analyse figures of merit after every n patterns\n"
"                             Default: 1000.  A value of zero means to do the\n"
"                             analysis only after reading all the patterns.\n"
"  -r, --rvsq                Output lists of R vs |q| (\"Luzzatti plots\")\n"
"                             when analysing figures of merit.\n"
"      --zone-axis           Output an [001] zone axis pattern each time the\n"
65
66
67
68
69
"                             figures of merit are analysed.\n"
"      --detwin              Correlate each new pattern with the current\n"
"                             model and choose the best fitting out of the\n"
"                             allowable twins.\n"
);
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
}


static void write_RvsQ(const char *name, double *ref, double *trueref,
                       unsigned int *counts, double scale, UnitCell *cell)
{
	FILE *fh;
	double smax, sbracket;
	signed int h, k, l;

	fh = fopen(name, "w");

	smax = 0.0;
	for ( h=-INDMAX; h<INDMAX; h++ ) {
	for ( k=-INDMAX; k<INDMAX; k++ ) {
	for ( l=-INDMAX; l<INDMAX; l++ ) {
Thomas White's avatar
Thomas White committed
86
		double s = 2.0*resolution(cell, h, k, l);
87
88
89
90
91
92
93
		if ( (lookup_count(counts, h, k, l) > 0) && (s > smax) ) {
			smax = s;
		}
	}
	}
	}

Thomas White's avatar
Thomas White committed
94
	for ( sbracket=0.0; sbracket<smax; sbracket+=smax/RVQDV ) {
95
96
97

		double top = 0.0;
		double bot = 0.0;
Thomas White's avatar
Thomas White committed
98
99
100
101
		int nhits = 0;
		int nrefl = 0;
		double R;
		double hits_per_refl;
102
103
104
105
106
107
108

		for ( h=-INDMAX; h<INDMAX; h++ ) {
		for ( k=-INDMAX; k<INDMAX; k++ ) {
		for ( l=-INDMAX; l<INDMAX; l++ ) {

			double s;
			int c;
Thomas White's avatar
Thomas White committed
109
110
111

			if ( (h==0) && (k==0) && (l==0) ) continue;

112
			c = lookup_count(counts, h, k, l);
Thomas White's avatar
Thomas White committed
113
			s = 2.0*resolution(cell, h, k, l);
Thomas White's avatar
Thomas White committed
114
			if ((s>=sbracket) && (s<sbracket+smax/RVQDV) && (c>0)) {
115
116
117
118
119
120
121

				double obs, calc, obsi;

				obs = lookup_intensity(ref, h, k, l);
				calc = lookup_intensity(trueref, h, k, l);

				obsi = obs / (double)c;
Thomas White's avatar
Thomas White committed
122
123
124
125
126
				top += pow(obsi - scale*calc, 2.0);
				bot += pow(obsi, 2.0);
				nhits += c;
				nrefl++;

127
128
129
130
131
132
			}

		}
		}
		}

Thomas White's avatar
Thomas White committed
133
134
135
136
137
		R = sqrt(top/bot);
		hits_per_refl = nrefl ? (double)nhits/nrefl : 0;

		fprintf(fh, "%8.5f %8.5f %5.2f\n", sbracket+smax/(2.0*RVQDV),
		                                R, hits_per_refl);
138
139
140
141
142
143
144
145

	}
	fclose(fh);
}


static void process_reflections(double *ref, double *trueref,
                                unsigned int *counts, unsigned int n_patterns,
Thomas White's avatar
Thomas White committed
146
                                UnitCell *cell, int do_rvsq, int do_zoneaxis)
147
148
149
150
151
152
{
	int j;
	double mean_counts;
	int ctot = 0;
	int nmeas = 0;
	double R, scale;
153
	FILE *fh;
154
155
156
157
158
159
160
161

	for ( j=0; j<LIST_SIZE; j++ ) {
		ctot += counts[j];
		if ( counts[j] > 0 ) nmeas++;
	}
	mean_counts = (double)ctot/nmeas;

	R = stat_r2(ref, trueref, counts, LIST_SIZE, &scale);
Thomas White's avatar
Thomas White committed
162
	STATUS("%8u: R=%5.2f%% (sf=%7.4e)  mean meas/refl=%5.2f,"
163
164
	       " %i reflections measured\n",
	       n_patterns, R*100.0, scale, mean_counts, nmeas);
165

Thomas White's avatar
Thomas White committed
166
167
	if ( do_rvsq ) {
		/* Record graph of R against q for this N */
Thomas White's avatar
Thomas White committed
168
		char name[64];
169
		snprintf(name, 63, "R_vs_q-%u.dat", n_patterns);
Thomas White's avatar
Thomas White committed
170
171
		write_RvsQ(name, ref, trueref, counts, scale, cell);
	}
Thomas White's avatar
Thomas White committed
172
173
174

	if ( do_zoneaxis ) {
		char name[64];
175
		snprintf(name, 63, "ZA-%u.dat", n_patterns);
176
		write_reflections(name, counts, ref, 1, cell, 1);
Thomas White's avatar
Thomas White committed
177
	}
178
179
180
181

	fh = fopen("results/convergence.dat", "a");
	fprintf(fh, "%u %5.2f %5.2f\n", n_patterns, R*100.0, mean_counts);
	fclose(fh);
182
183
184
}


185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
static void merge_pattern(double *model, const double *new,
                          unsigned int *model_counts,
                          const unsigned int *counts, int mo, int sum)
{
	signed int h, k, l;

	for ( l=-INDMAX; l<INDMAX; l++ ) {
	for ( k=-INDMAX; k<INDMAX; k++ ) {
	for ( h=-INDMAX; h<INDMAX; h++ ) {

		double intensity;

		if ( lookup_count(counts, h, k, l) == 0 ) continue;

		intensity = lookup_intensity(new, h, k, l);

		if ( !mo ) {
			integrate_intensity(model, h, k, l, intensity);
			if ( sum ) {
				set_count(model_counts, h, k, l, 1);
			} else {
				integrate_count(model_counts, h, k, l, 1);
			}
		} else {
			if ( intensity > lookup_intensity(model, h, k, l) ) {
				set_intensity(model, h, k, l, intensity);
			}
			set_count(model_counts, h, k, l, 1);
		}

	}
	}
	}
}


221
222
223
224
int main(int argc, char *argv[])
{
	int c;
	char *filename = NULL;
225
	char *output = NULL;
226
227
	FILE *fh;
	unsigned int n_patterns;
228
229
	double *model, *trueref = NULL;
	unsigned int *model_counts;
230
	char *rval;
231
	UnitCell *cell;
Thomas White's avatar
Thomas White committed
232
	int config_maxonly = 0;
Thomas White's avatar
Thomas White committed
233
	int config_every = 1000;
Thomas White's avatar
Thomas White committed
234
	int config_rvsq = 0;
Thomas White's avatar
Thomas White committed
235
	int config_stopafter = 0;
Thomas White's avatar
Thomas White committed
236
	int config_zoneaxis = 0;
237
	int config_sum = 0;
238
	int config_detwin = 0;
239
	char *intfile = NULL;
240
241
242
	double *new_pattern = NULL;
	unsigned int *new_counts = NULL;
	unsigned int n_total_patterns;
243
244
245
246
247

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
248
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
249
		{"max-only",           0, &config_maxonly,     1},
Thomas White's avatar
Thomas White committed
250
		{"output-every",       1, NULL,               'e'},
Thomas White's avatar
Thomas White committed
251
		{"rvsq",               0, NULL,               'r'},
Thomas White's avatar
Thomas White committed
252
		{"stop-after",         1, NULL,               's'},
Thomas White's avatar
Thomas White committed
253
		{"zone-axis",          0, &config_zoneaxis,    1},
254
		{"compare-with",       0, NULL,               'c'},
255
256
		{"sum",                0, &config_sum,         1},
		{"detwin",             0, &config_detwin,      1},
257
258
259
260
		{0, 0, NULL, 0}
	};

	/* Short options */
261
	while ((c = getopt_long(argc, argv, "hi:e:ro:", longopts, NULL)) != -1) {
262
263
264
265
266
267
268
269
270
271
272
273

		switch (c) {
		case 'h' : {
			show_help(argv[0]);
			return 0;
		}

		case 'i' : {
			filename = strdup(optarg);
			break;
		}

274
275
276
277
278
		case 'o' : {
			output = strdup(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
279
280
281
282
283
		case 'r' : {
			config_rvsq = 1;
			break;
		}

Thomas White's avatar
Thomas White committed
284
285
286
287
288
		case 'e' : {
			config_every = atoi(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
289
290
291
292
293
		case 's' : {
			config_stopafter = atoi(optarg);
			break;
		}

294
295
296
297
298
		case 'c' : {
			intfile = strdup(optarg);
			break;
		}

299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
		case 0 : {
			break;
		}

		default : {
			return 1;
		}
		}

	}

	if ( filename == NULL ) {
		ERROR("Please specify filename using the -i option\n");
		return 1;
	}

315
316
	if ( intfile != NULL ) {
		STATUS("Comparing against '%s'\n", intfile);
317
		trueref = read_reflections(intfile, NULL);
318
319
320
321
322
		free(intfile);
	} else {
		trueref = NULL;
	}

323
324
	model = new_list_intensity();
	model_counts = new_list_count();
325
	cell = load_cell_from_pdb("molecule.pdb");
326
327
	new_pattern = new_list_intensity();
	new_counts = new_list_count();
328
329
330
331
332
333
334
335
336
337
338
339

	if ( strcmp(filename, "-") == 0 ) {
		fh = stdin;
	} else {
		fh = fopen(filename, "r");
	}
	free(filename);
	if ( fh == NULL ) {
		ERROR("Failed to open input file\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
340
341
342
343
344
345
346
347
348
349
350
351
352
353
	/* Count the number of patterns in the file */
	n_total_patterns = 0;
	do {
		char line[1024];

		rval = fgets(line, 1023, fh);
		if ( (strncmp(line, "Reflections from indexing", 25) == 0)
		    || (strncmp(line, "New pattern", 11) == 0) ) {
		    n_total_patterns++;
		}
	} while ( rval != NULL );
	rewind(fh);
	STATUS("There are %i patterns to process\n", n_total_patterns);

354
355
356
357
	n_patterns = 0;
	do {

		char line[1024];
358
359
		signed int h, k, l;
		float intensity;
360
361
362
		int r;

		rval = fgets(line, 1023, fh);
363
364
		if ( (strncmp(line, "Reflections from indexing", 25) == 0)
		    || (strncmp(line, "New pattern", 11) == 0) ) {
Thomas White's avatar
Tidy-up    
Thomas White committed
365

366
			/* Start of first pattern? */
Thomas White's avatar
Thomas White committed
367
368
369
370
			if ( n_patterns == 0 ) {
				n_patterns++;
				continue;
			}
Thomas White's avatar
Tidy-up    
Thomas White committed
371

372
373
374
375
376
			if ( config_detwin ) {
				detwin_intensities(model, new_pattern,
				                   model_counts, new_counts);
			}

377
378
379
380
			/* Start of second or later pattern */
			merge_pattern(model, new_pattern, model_counts,
			              new_counts, config_maxonly, config_sum);

381
			if (config_every && (n_patterns % config_every == 0)) {
382
383
384
				process_reflections(model, trueref,
				                    model_counts, n_patterns,
				                    cell, config_rvsq,
Thomas White's avatar
Thomas White committed
385
				                    config_zoneaxis);
386
			}
Thomas White's avatar
Thomas White committed
387
388
389

			if ( n_patterns == config_stopafter ) break;

390
391
			zero_list_count(new_counts);

Thomas White's avatar
Thomas White committed
392
			n_patterns++;
Thomas White's avatar
Thomas White committed
393
394
395

			progress_bar(n_patterns, n_total_patterns, "Merging");

396
397
		}

398
		r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity);
399
400
		if ( r != 4 ) continue;

Thomas White's avatar
Tidy-up    
Thomas White committed
401
402
		if ( (h==0) && (k==0) && (l==0) ) continue;

403
404
405
		if ( lookup_count(new_counts, h, k, l) != 0 ) {
			ERROR("More than one measurement for %i %i %i in"
			      " pattern number %i\n", h, k, l, n_patterns);
Thomas White's avatar
Thomas White committed
406
		}
407
408
		set_intensity(new_pattern, h, k, l, intensity);
		set_count(new_counts, h, k, l, 1);
409
410
411
412
413

	} while ( rval != NULL );

	fclose(fh);

414
	if ( trueref != NULL ) {
415
416
		process_reflections(model, trueref, model_counts, n_patterns,
		                    cell, config_rvsq, config_zoneaxis);
Thomas White's avatar
Thomas White committed
417
	}
418

419
	if ( output != NULL ) {
420
		write_reflections(output, model_counts, model, 0, cell, 1);
421
	}
Thomas White's avatar
Thomas White committed
422

423
424
425
	if ( config_zoneaxis ) {
		char name[64];
		snprintf(name, 63, "ZA-%u.dat", n_patterns);
426
		write_reflections(name, model_counts, model, 1, cell, 10);
427
428
	}

Thomas White's avatar
Thomas White committed
429
	STATUS("There were %u patterns.\n", n_patterns);
430
431
432

	return 0;
}