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
static void merge_pattern(double *model, const double *new,
                          unsigned int *model_counts,
191
                          ReflItemList *items, int mo, int sum)
192
{
193
	int i;
194

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

		double intensity;
198
199
		signed int h, k, l;
		struct refl_item *item;
200

201
202
203
204
205
		item = get_item(items, i);

		h = item->h;
		k = item->k;
		l = item->l;
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226

		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);
		}

	}
}


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

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

	/* Short options */
277
278
	while ((c = getopt_long(argc, argv, "hi:e:ro:p:y:",
	                        longopts, NULL)) != -1) {
279
280

		switch (c) {
Thomas White's avatar
Thomas White committed
281
		case 'h' :
282
283
284
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
285
		case 'i' :
286
287
288
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
289
		case 'o' :
290
291
292
			output = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
293
294

		case 'r' :
Thomas White's avatar
Thomas White committed
295
296
297
			config_rvsq = 1;
			break;

Thomas White's avatar
Thomas White committed
298
		case 'e' :
Thomas White's avatar
Thomas White committed
299
300
301
			config_every = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
302
		case 's' :
Thomas White's avatar
Thomas White committed
303
304
305
			config_stopafter = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
306
		case 'c' :
307
308
309
			intfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
310
		case 'p' :
311
312
313
			pdb = strdup(optarg);
			break;

314
315
316
317
		case 'y' :
			sym = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
318
		case 0 :
319
320
			break;

Thomas White's avatar
Thomas White committed
321
		default :
322
323
324
325
326
327
328
329
330
331
			return 1;
		}

	}

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

332
	if ( intfile != NULL ) {
333
		truecounts = new_list_count();
334
		STATUS("Comparing against '%s'\n", intfile);
335
		trueref = read_reflections(intfile, truecounts, NULL);
336
337
338
339
340
		free(intfile);
	} else {
		trueref = NULL;
	}

341
342
343
344
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

345
346
	model = new_list_intensity();
	model_counts = new_list_count();
347
	cell = load_cell_from_pdb(pdb);
348
	new_pattern = new_list_intensity();
349
	items = new_items();
350
351
352
353
354
355
356
357
358
359
360
361

	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
362
363
364
365
366
367
368
369
370
371
372
373
374
375
	/* 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);

376
	n_patterns = 0;
Thomas White's avatar
Thomas White committed
377
	f0_valid = 0;
378
379
380
	do {

		char line[1024];
381
382
		signed int h, k, l;
		float intensity;
383
384
385
		int r;

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

389
			/* Start of first pattern? */
Thomas White's avatar
Thomas White committed
390
391
392
393
			if ( n_patterns == 0 ) {
				n_patterns++;
				continue;
			}
Thomas White's avatar
Tidy-up    
Thomas White committed
394

395
			/* Detwin before scaling */
396
397
			if ( config_detwin ) {
				detwin_intensities(model, new_pattern,
398
				                   model_counts, items);
399
400
			}

Thomas White's avatar
Thomas White committed
401
402
			/* Assume a default I0 if we don't have one by now */
			if ( config_scale && !f0_valid ) {
403
				n_nof0++;
404
405
406
				f0 = 1.0;
			}

407
408
409
			/* Scale if requested */
			if ( config_scale ) {
				scale_intensities(model, new_pattern,
410
				                  model_counts, items, f0,
Thomas White's avatar
Thomas White committed
411
				                  f0_valid);
412
413
			}

414
415
			/* Start of second or later pattern */
			merge_pattern(model, new_pattern, model_counts,
416
			              items, config_maxonly, config_sum);
417

418
419
			if ( (trueref != NULL) && config_every
			    && (n_patterns % config_every == 0) ) {
420
421
422
423
				process_reflections(model, model_counts,
				                    trueref, truecounts,
				                    n_patterns, cell,
				                    config_rvsq,
Thomas White's avatar
Thomas White committed
424
				                    config_zoneaxis);
425
			}
Thomas White's avatar
Thomas White committed
426
427
428

			if ( n_patterns == config_stopafter ) break;

Thomas White's avatar
Thomas White committed
429
			n_patterns++;
430
			clear_items(items);
Thomas White's avatar
Thomas White committed
431
432
433

			progress_bar(n_patterns, n_total_patterns, "Merging");

Thomas White's avatar
Thomas White committed
434
			f0_valid = 0;
435
436
437
438
439
440

		}

		if ( strncmp(line, "f0 = ", 5) == 0 ) {
			r = sscanf(line, "f0 = %f", &f0);
			if ( r != 1 ) {
Thomas White's avatar
Thomas White committed
441
442
				f0 = 1.0;
				f0_valid = 0;
443
444
				continue;
			}
Thomas White's avatar
Thomas White committed
445
			f0_valid = 1;
446
447
		}

448
		r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity);
449
450
		if ( r != 4 ) continue;

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

453
		if ( find_item(items, h, k, l) != 0 ) {
454
455
			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
456
		}
457
		set_intensity(new_pattern, h, k, l, intensity);
458
		add_item(items, h, k, l);
459
460
461
462
463

	} while ( rval != NULL );

	fclose(fh);

464
	if ( trueref != NULL ) {
465
466
467
		process_reflections(model, model_counts, trueref, truecounts,
		                    n_patterns, cell, config_rvsq,
		                    config_zoneaxis);
Thomas White's avatar
Thomas White committed
468
	}
469

470
	if ( output != NULL ) {
471
472
		write_reflections(output, model_counts, model, NULL,
		                  0, cell, 1);
473
	}
Thomas White's avatar
Thomas White committed
474

475
476
477
	if ( config_zoneaxis ) {
		char name[64];
		snprintf(name, 63, "ZA-%u.dat", n_patterns);
478
		write_reflections(name, model_counts, model, NULL, 1, cell, 10);
479
480
	}

Thomas White's avatar
Thomas White committed
481
	STATUS("There were %u patterns.\n", n_patterns);
482
	STATUS("%i had no f0 valid value.\n", n_nof0);
483

484
485
	delete_items(items);

486
487
	return 0;
}