process_hkl.c 9.96 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
29


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


34
35
36
static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
Thomas White's avatar
Thomas White committed
37
38
39
	printf(
"Assemble and process FEL Bragg intensities.\n"
"\n"
40
41
42
43
"  -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
44
"\n"
45
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"
"                             figures of merit are analysed.\n");
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
}


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
81
		double s = 2.0*resolution(cell, h, k, l);
82
83
84
85
86
87
88
		if ( (lookup_count(counts, h, k, l) > 0) && (s > smax) ) {
			smax = s;
		}
	}
	}
	}

Thomas White's avatar
Thomas White committed
89
	for ( sbracket=0.0; sbracket<smax; sbracket+=smax/RVQDV ) {
90
91
92

		double top = 0.0;
		double bot = 0.0;
Thomas White's avatar
Thomas White committed
93
94
95
96
		int nhits = 0;
		int nrefl = 0;
		double R;
		double hits_per_refl;
97
98
99
100
101
102
103

		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
104
105
106

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

107
			c = lookup_count(counts, h, k, l);
Thomas White's avatar
Thomas White committed
108
			s = 2.0*resolution(cell, h, k, l);
Thomas White's avatar
Thomas White committed
109
			if ((s>=sbracket) && (s<sbracket+smax/RVQDV) && (c>0)) {
110
111
112
113
114
115
116

				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
117
118
119
120
121
				top += pow(obsi - scale*calc, 2.0);
				bot += pow(obsi, 2.0);
				nhits += c;
				nrefl++;

122
123
124
125
126
127
			}

		}
		}
		}

Thomas White's avatar
Thomas White committed
128
129
130
131
132
		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);
133
134
135
136
137
138
139
140

	}
	fclose(fh);
}


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

	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
157
	STATUS("%8u: R=%5.2f%% (sf=%7.4e)  mean meas/refl=%5.2f,"
158
159
	       " %i reflections measured\n",
	       n_patterns, R*100.0, scale, mean_counts, nmeas);
160

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

	if ( do_zoneaxis ) {
		char name[64];
170
		snprintf(name, 63, "ZA-%u.dat", n_patterns);
171
		write_reflections(name, counts, ref, 1, cell, 1);
Thomas White's avatar
Thomas White committed
172
	}
173
174
175
176

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


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

	}
	}
	}
}


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

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
242
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
243
		{"max-only",           0, &config_maxonly,     1},
Thomas White's avatar
Thomas White committed
244
		{"output-every",       1, NULL,               'e'},
Thomas White's avatar
Thomas White committed
245
		{"rvsq",               0, NULL,               'r'},
Thomas White's avatar
Thomas White committed
246
		{"stop-after",         1, NULL,               's'},
Thomas White's avatar
Thomas White committed
247
		{"zone-axis",          0, &config_zoneaxis,    1},
248
		{"compare-with",       0, NULL,               'c'},
249
		{"sum",                0, &config_sum,    1},
250
251
252
253
		{0, 0, NULL, 0}
	};

	/* Short options */
254
	while ((c = getopt_long(argc, argv, "hi:e:ro:", longopts, NULL)) != -1) {
255
256
257
258
259
260
261
262
263
264
265
266

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

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

267
268
269
270
271
		case 'o' : {
			output = strdup(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
272
273
274
275
276
		case 'r' : {
			config_rvsq = 1;
			break;
		}

Thomas White's avatar
Thomas White committed
277
278
279
280
281
		case 'e' : {
			config_every = atoi(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
282
283
284
285
286
		case 's' : {
			config_stopafter = atoi(optarg);
			break;
		}

287
288
289
290
291
		case 'c' : {
			intfile = strdup(optarg);
			break;
		}

292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
		case 0 : {
			break;
		}

		default : {
			return 1;
		}
		}

	}

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

308
309
	if ( intfile != NULL ) {
		STATUS("Comparing against '%s'\n", intfile);
310
		trueref = read_reflections(intfile, NULL);
311
312
313
314
315
		free(intfile);
	} else {
		trueref = NULL;
	}

316
317
	model = new_list_intensity();
	model_counts = new_list_count();
318
	cell = load_cell_from_pdb("molecule.pdb");
319
320
	new_pattern = new_list_intensity();
	new_counts = new_list_count();
321
322
323
324
325
326
327
328
329
330
331
332

	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
333
334
335
336
337
338
339
340
341
342
343
344
345
346
	/* 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);

347
348
349
350
	n_patterns = 0;
	do {

		char line[1024];
351
352
		signed int h, k, l;
		float intensity;
353
354
355
		int r;

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

359
			/* Start of first pattern? */
Thomas White's avatar
Thomas White committed
360
361
362
363
			if ( n_patterns == 0 ) {
				n_patterns++;
				continue;
			}
Thomas White's avatar
Tidy-up    
Thomas White committed
364

365
366
367
368
			/* Start of second or later pattern */
			merge_pattern(model, new_pattern, model_counts,
			              new_counts, config_maxonly, config_sum);

369
			if (config_every && (n_patterns % config_every == 0)) {
370
371
372
				process_reflections(model, trueref,
				                    model_counts, n_patterns,
				                    cell, config_rvsq,
Thomas White's avatar
Thomas White committed
373
				                    config_zoneaxis);
374
			}
Thomas White's avatar
Thomas White committed
375
376
377

			if ( n_patterns == config_stopafter ) break;

378
379
			zero_list_count(new_counts);

Thomas White's avatar
Thomas White committed
380
			n_patterns++;
Thomas White's avatar
Thomas White committed
381
382
383

			progress_bar(n_patterns, n_total_patterns, "Merging");

384
385
		}

386
		r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity);
387
388
		if ( r != 4 ) continue;

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

391
392
393
		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
394
		}
395
396
		set_intensity(new_pattern, h, k, l, intensity);
		set_count(new_counts, h, k, l, 1);
397
398
399
400
401

	} while ( rval != NULL );

	fclose(fh);

402
	if ( trueref != NULL ) {
403
404
		process_reflections(model, trueref, model_counts, n_patterns,
		                    cell, config_rvsq, config_zoneaxis);
Thomas White's avatar
Thomas White committed
405
	}
406

407
	if ( output != NULL ) {
408
		write_reflections(output, model_counts, model, 0, cell, 1);
409
	}
Thomas White's avatar
Thomas White committed
410

411
412
413
	if ( config_zoneaxis ) {
		char name[64];
		snprintf(name, 63, "ZA-%u.dat", n_patterns);
414
		write_reflections(name, model_counts, model, 1, cell, 10);
415
416
	}

Thomas White's avatar
Thomas White committed
417
	STATUS("There were %u patterns.\n", n_patterns);
418
419
420

	return 0;
}