process_hkl.c 11.8 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 *sym = NULL;
250
	char *pdb = NULL;
251
	float f0;
Thomas White's avatar
Thomas White committed
252
	int f0_valid;
253
	int n_nof0 = 0;
254
255
256
257
258

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

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

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

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

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

Thomas White's avatar
Thomas White committed
291
292

		case 'r' :
Thomas White's avatar
Thomas White committed
293
294
295
			config_rvsq = 1;
			break;

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

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

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

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

312
313
314
315
		case 'y' :
			sym = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
316
		case 0 :
317
318
			break;

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

	}

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

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

339
340
341
342
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

343
344
	model = new_list_intensity();
	model_counts = new_list_count();
345
	cell = load_cell_from_pdb(pdb);
346
347
	new_pattern = new_list_intensity();
	new_counts = new_list_count();
348
349
350
351
352
353
354
355
356
357
358
359

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

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

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

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

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

393
			/* Detwin before scaling */
394
395
396
397
398
			if ( config_detwin ) {
				detwin_intensities(model, new_pattern,
				                   model_counts, new_counts);
			}

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

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

412
413
414
415
			/* Start of second or later pattern */
			merge_pattern(model, new_pattern, model_counts,
			              new_counts, config_maxonly, config_sum);

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

			if ( n_patterns == config_stopafter ) break;

427
428
			zero_list_count(new_counts);

Thomas White's avatar
Thomas White committed
429
			n_patterns++;
Thomas White's avatar
Thomas White committed
430
431
432

			progress_bar(n_patterns, n_total_patterns, "Merging");

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

		}

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

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

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

452
453
454
		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
455
		}
456
457
		set_intensity(new_pattern, h, k, l, intensity);
		set_count(new_counts, h, k, l, 1);
458
459
460
461
462

	} while ( rval != NULL );

	fclose(fh);

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

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

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

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

	return 0;
}