process_hkl.c 11.5 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
"                             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"
69
70
"      --scale               Scale each pattern for best fit with the current\n"
"                             model.\n"
71
);
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
}


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

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

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

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

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

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

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

129
130
131
132
133
134
			}

		}
		}
		}

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

	}
	fclose(fh);
}


146
147
148
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
149
                                UnitCell *cell, int do_rvsq, int do_zoneaxis)
150
151
152
153
154
155
{
	int j;
	double mean_counts;
	int ctot = 0;
	int nmeas = 0;
	double R, scale;
156
	FILE *fh;
157
158
159
160
161
162
163

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

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

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

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

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

	}
	}
	}
}


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

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

	/* Short options */
269
	while ((c = getopt_long(argc, argv, "hi:e:ro:", longopts, NULL)) != -1) {
270
271
272
273
274
275
276
277
278
279
280
281

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

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

282
283
284
285
286
		case 'o' : {
			output = strdup(optarg);
			break;
		}

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

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

Thomas White's avatar
Thomas White committed
297
298
299
300
301
		case 's' : {
			config_stopafter = atoi(optarg);
			break;
		}

302
303
304
305
306
		case 'c' : {
			intfile = strdup(optarg);
			break;
		}

307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
		case 0 : {
			break;
		}

		default : {
			return 1;
		}
		}

	}

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

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

332
333
	model = new_list_intensity();
	model_counts = new_list_count();
334
	cell = load_cell_from_pdb("molecule.pdb");
335
336
	new_pattern = new_list_intensity();
	new_counts = new_list_count();
337
338
339
340
341
342
343
344
345
346
347
348

	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
349
350
351
352
353
354
355
356
357
358
359
360
361
362
	/* 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);

363
	n_patterns = 0;
Thomas White's avatar
Thomas White committed
364
	f0_valid = 0;
365
366
367
	do {

		char line[1024];
368
369
		signed int h, k, l;
		float intensity;
370
371
372
		int r;

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

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

382
			/* Detwin before scaling */
383
384
385
386
387
			if ( config_detwin ) {
				detwin_intensities(model, new_pattern,
				                   model_counts, new_counts);
			}

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

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

401
402
403
404
			/* Start of second or later pattern */
			merge_pattern(model, new_pattern, model_counts,
			              new_counts, config_maxonly, config_sum);

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

			if ( n_patterns == config_stopafter ) break;

416
417
			zero_list_count(new_counts);

Thomas White's avatar
Thomas White committed
418
			n_patterns++;
Thomas White's avatar
Thomas White committed
419
420
421

			progress_bar(n_patterns, n_total_patterns, "Merging");

Thomas White's avatar
Thomas White committed
422
			f0_valid = 0;
423
424
425
426
427
428
429

		}

		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
430
431
				f0 = 1.0;
				f0_valid = 0;
432
433
				continue;
			}
Thomas White's avatar
Thomas White committed
434
			f0_valid = 1;
435
436
		}

437
		r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity);
438
439
		if ( r != 4 ) continue;

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

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

	} while ( rval != NULL );

	fclose(fh);

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

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

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

Thomas White's avatar
Thomas White committed
470
	STATUS("There were %u patterns.\n", n_patterns);
471
472
473

	return 0;
}