process_hkl.c 11.7 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"
45
"  -p, --pdb=<filename>      PDB file to use (default: molecule.pdb).\n"
Thomas White's avatar
Thomas White committed
46
"\n"
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
"      --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"
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
"      --scale               Scale each pattern for best fit with the current\n"
"                             model.\n"
72
);
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
}


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
89
		double s = 2.0*resolution(cell, h, k, l);
90
91
92
93
94
95
96
		if ( (lookup_count(counts, h, k, l) > 0) && (s > smax) ) {
			smax = s;
		}
	}
	}
	}

Thomas White's avatar
Thomas White committed
97
	for ( sbracket=0.0; sbracket<smax; sbracket+=smax/RVQDV ) {
98
99
100

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

		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
112
113
114

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

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

				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
125
126
127
128
129
				top += pow(obsi - scale*calc, 2.0);
				bot += pow(obsi, 2.0);
				nhits += c;
				nrefl++;

130
131
132
133
134
135
			}

		}
		}
		}

Thomas White's avatar
Thomas White committed
136
137
138
139
140
		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);
141
142
143
144
145
146

	}
	fclose(fh);
}


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

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

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

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

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

	fh = fopen("results/convergence.dat", "a");
	fprintf(fh, "%u %5.2f %5.2f\n", n_patterns, R*100.0, mean_counts);
	fclose(fh);
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
221
222
223
224
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);
		}

	}
	}
	}
}


225
226
227
228
int main(int argc, char *argv[])
{
	int c;
	char *filename = NULL;
229
	char *output = NULL;
230
231
	FILE *fh;
	unsigned int n_patterns;
232
233
	double *model, *trueref = NULL;
	unsigned int *model_counts;
234
	char *rval;
235
	UnitCell *cell;
Thomas White's avatar
Thomas White committed
236
	int config_maxonly = 0;
Thomas White's avatar
Thomas White committed
237
	int config_every = 1000;
Thomas White's avatar
Thomas White committed
238
	int config_rvsq = 0;
Thomas White's avatar
Thomas White committed
239
	int config_stopafter = 0;
Thomas White's avatar
Thomas White committed
240
	int config_zoneaxis = 0;
241
	int config_sum = 0;
242
	int config_detwin = 0;
243
	int config_scale = 0;
244
	char *intfile = NULL;
245
246
247
	double *new_pattern = NULL;
	unsigned int *new_counts = NULL;
	unsigned int n_total_patterns;
248
	unsigned int *truecounts = NULL;
249
	char *pdb = NULL;
250
	float f0;
Thomas White's avatar
Thomas White committed
251
	int f0_valid;
252
253
254
255
256

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
257
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
258
		{"max-only",           0, &config_maxonly,     1},
Thomas White's avatar
Thomas White committed
259
		{"output-every",       1, NULL,               'e'},
Thomas White's avatar
Thomas White committed
260
		{"rvsq",               0, NULL,               'r'},
Thomas White's avatar
Thomas White committed
261
		{"stop-after",         1, NULL,               's'},
Thomas White's avatar
Thomas White committed
262
		{"zone-axis",          0, &config_zoneaxis,    1},
263
		{"compare-with",       0, NULL,               'c'},
264
265
		{"sum",                0, &config_sum,         1},
		{"detwin",             0, &config_detwin,      1},
266
267
		{"scale",              0, &config_scale,       1},
		{"pdb",                1, NULL,               'p'},
268
269
270
271
		{0, 0, NULL, 0}
	};

	/* Short options */
272
	while ((c = getopt_long(argc, argv, "hi:e:ro:p:", longopts, NULL)) != -1) {
273
274

		switch (c) {
Thomas White's avatar
Thomas White committed
275
		case 'h' :
276
277
278
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
279
		case 'i' :
280
281
282
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
283
		case 'o' :
284
285
286
			output = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
287
288

		case 'r' :
Thomas White's avatar
Thomas White committed
289
290
291
			config_rvsq = 1;
			break;

Thomas White's avatar
Thomas White committed
292
		case 'e' :
Thomas White's avatar
Thomas White committed
293
294
295
			config_every = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
296
		case 's' :
Thomas White's avatar
Thomas White committed
297
298
299
			config_stopafter = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
300
		case 'c' :
301
302
303
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
304
		case 'p' :
305
306
307
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
308
		case 0 :
309
310
			break;

Thomas White's avatar
Thomas White committed
311
		default :
312
313
314
315
316
317
318
319
320
321
			return 1;
		}

	}

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

322
	if ( intfile != NULL ) {
323
		truecounts = new_list_count();
324
		STATUS("Comparing against '%s'\n", intfile);
325
		trueref = read_reflections(intfile, truecounts, NULL);
326
327
328
329
330
		free(intfile);
	} else {
		trueref = NULL;
	}

331
332
333
334
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

335
336
	model = new_list_intensity();
	model_counts = new_list_count();
337
	cell = load_cell_from_pdb(pdb);
338
339
	new_pattern = new_list_intensity();
	new_counts = new_list_count();
340
341
342
343
344
345
346
347
348
349
350
351

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

366
	n_patterns = 0;
Thomas White's avatar
Thomas White committed
367
	f0_valid = 0;
368
369
370
	do {

		char line[1024];
371
372
		signed int h, k, l;
		float intensity;
373
374
375
		int r;

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

379
			/* Start of first pattern? */
Thomas White's avatar
Thomas White committed
380
381
382
383
			if ( n_patterns == 0 ) {
				n_patterns++;
				continue;
			}
Thomas White's avatar
Tidy-up    
Thomas White committed
384

385
			/* Detwin before scaling */
386
387
388
389
390
			if ( config_detwin ) {
				detwin_intensities(model, new_pattern,
				                   model_counts, new_counts);
			}

Thomas White's avatar
Thomas White committed
391
392
			/* Assume a default I0 if we don't have one by now */
			if ( config_scale && !f0_valid ) {
393
394
395
396
				ERROR("No f0 value.\n");
				f0 = 1.0;
			}

397
398
399
			/* Scale if requested */
			if ( config_scale ) {
				scale_intensities(model, new_pattern,
Thomas White's avatar
Thomas White committed
400
401
				                  model_counts, new_counts, f0,
				                  f0_valid);
402
403
			}

404
405
406
407
			/* Start of second or later pattern */
			merge_pattern(model, new_pattern, model_counts,
			              new_counts, config_maxonly, config_sum);

408
409
			if ( (trueref != NULL) && config_every
			    && (n_patterns % config_every == 0) ) {
410
411
412
413
				process_reflections(model, model_counts,
				                    trueref, truecounts,
				                    n_patterns, cell,
				                    config_rvsq,
Thomas White's avatar
Thomas White committed
414
				                    config_zoneaxis);
415
			}
Thomas White's avatar
Thomas White committed
416
417
418

			if ( n_patterns == config_stopafter ) break;

419
420
			zero_list_count(new_counts);

Thomas White's avatar
Thomas White committed
421
			n_patterns++;
Thomas White's avatar
Thomas White committed
422
423
424

			progress_bar(n_patterns, n_total_patterns, "Merging");

Thomas White's avatar
Thomas White committed
425
			f0_valid = 0;
426
427
428
429
430
431
432

		}

		if ( strncmp(line, "f0 = ", 5) == 0 ) {
			r = sscanf(line, "f0 = %f", &f0);
			if ( r != 1 ) {
				ERROR("Couldn't understand f0 line.\n");
Thomas White's avatar
Thomas White committed
433
434
				f0 = 1.0;
				f0_valid = 0;
435
436
				continue;
			}
Thomas White's avatar
Thomas White committed
437
			f0_valid = 1;
438
439
		}

440
		r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity);
441
442
		if ( r != 4 ) continue;

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

445
446
447
		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
448
		}
449
450
		set_intensity(new_pattern, h, k, l, intensity);
		set_count(new_counts, h, k, l, 1);
451
452
453
454
455

	} while ( rval != NULL );

	fclose(fh);

456
	if ( trueref != NULL ) {
457
458
459
		process_reflections(model, model_counts, trueref, truecounts,
		                    n_patterns, cell, config_rvsq,
		                    config_zoneaxis);
Thomas White's avatar
Thomas White committed
460
	}
461

462
	if ( output != NULL ) {
463
464
		write_reflections(output, model_counts, model, NULL,
		                  0, cell, 1);
465
	}
Thomas White's avatar
Thomas White committed
466

467
468
469
	if ( config_zoneaxis ) {
		char name[64];
		snprintf(name, 63, "ZA-%u.dat", n_patterns);
470
		write_reflections(name, model_counts, model, NULL, 1, cell, 10);
471
472
	}

Thomas White's avatar
Thomas White committed
473
	STATUS("There were %u patterns.\n", n_patterns);
474
475
476

	return 0;
}