process_hkl.c 14.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
7
8
 * Copyright © 2012 Thomas White <taw@physics.org>
 * Copyright © 2012 Andrew Martin <andrew.martin@desy.de>
 * Copyright © 2012 Lorenzo Galli <lorenzo.galli@desy.de>
9
 *
Thomas White's avatar
Thomas White committed
10
11
12
13
14
15
16
17
18
19
20
21
22
23
 * This file is part of CrystFEL.
 *
 * CrystFEL is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * CrystFEL is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
 *
 */


#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"
41
#include "reflist-utils.h"
42
#include "symmetry.h"
43
#include "stream.h"
44
45
#include "reflist.h"
#include "image.h"
46
47
48
49
50


static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
Thomas White's avatar
Thomas White committed
51
52
53
	printf(
"Assemble and process FEL Bragg intensities.\n"
"\n"
54
55
56
"  -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"
57
"                             Default: processed.hkl).\n"
Thomas White's avatar
Thomas White committed
58
"  -y, --symmetry=<sym>      Merge according to point group <sym>.\n"
Thomas White's avatar
Thomas White committed
59
"\n"
60
61
62
63
64
65
66
67
"      --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"
68
"      --start-after=<n>     Skip n patterns at the start of the stream.\n"
Thomas White's avatar
Thomas White committed
69
"      --stop-after=<n>      Stop after processing n patterns.\n"
70
71
"  -g, --histogram=<h,k,l>   Calculate the histogram of measurements for this\n"
"                             reflection.\n"
72
"  -z, --hist-parameters     Set the range for the histogram and the number of\n"
Thomas White's avatar
Thomas White committed
73
"          =<min,max,nbins>   bins. \n"
74
"\n"
75
76
"      --scale               Scale each pattern for best fit with the current\n"
"                             model.\n"
77
"      --reference=<file>    Compare against intensities from <file> when\n"
Thomas White's avatar
Thomas White committed
78
"                             scaling. \n"
79
);
80
81
82
}


Thomas White's avatar
Thomas White committed
83
84
static void plot_histogram(double *vals, int n, float hist_min, float hist_max,
                           int nbins)
85
86
87
88
89
{
	int i;
	double max = -INFINITY;
	double min = +INFINITY;
	double step;
90
	int histo[nbins];
91
92
93
94
95
96
97
98
	FILE *fh;

	fh = fopen("histogram.dat", "w");
	if ( fh == NULL ) {
		ERROR("Couldn't open 'histogram.dat'\n");
		return;
	}

99
100
101
102
103
104
105
106
	if ( hist_min == hist_max ) {
		for ( i=0; i<n; i++ ) {
			if ( vals[i] > max ) max = vals[i];
			if ( vals[i] < min ) min = vals[i];
		}
	} else {
		min = hist_min;
		max = hist_max;
107
	}
Thomas White's avatar
Thomas White committed
108
	STATUS("min max nbins: %f %f %i\n", min, max, nbins);
109
	min--;  max++;
110

111
	for ( i=0; i<nbins; i++ ) {
112
113
114
		histo[i] = 0;
	}

115
	step = (max-min)/nbins;
116
117
118

	for ( i=0; i<n; i++ ) {
		int bin;
119
120
121
122
		if ( (vals[i] > min) && (vals[i] < max) ) {
			bin = (vals[i]-min)/step;
			histo[bin]++;
		}
123
124
	}

125
	for ( i=0; i<nbins; i++ ) {
126
127
128
129
130
131
132
		fprintf(fh, "%f %i\n", min+step*i, histo[i]);
	}

	fclose(fh);
}


133
static void merge_pattern(RefList *model, RefList *new, int max_only,
Thomas White's avatar
Thomas White committed
134
                          const SymOpList *sym,
135
                          double *hist_vals, signed int hist_h,
136
137
                          signed int hist_k, signed int hist_l, int *hist_n,
                          int pass)
138
{
139
140
	Reflection *refl;
	RefListIterator *iter;
141

142
143
144
	for ( refl = first_refl(new, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
145
146

		double intensity;
147
		signed int h, k, l;
148
149
		Reflection *model_version;
		double model_int;
150

151
		get_indices(refl, &h, &k, &l);
Thomas White's avatar
Thomas White committed
152
153

		/* Put into the asymmetric unit for the target group */
Thomas White's avatar
Thomas White committed
154
		get_asymm(sym, h, k, l, &h, &k, &l);
Thomas White's avatar
Thomas White committed
155

156
157
158
159
		model_version = find_refl(model, h, k, l);
		if ( model_version == NULL ) {
			model_version = add_refl(model, h, k, l);
		}
160

161
162
		/* Read the intensity from the original location
		 * (i.e. before screwing around with symmetry) */
163
164
165
166
		intensity = get_intensity(refl);

		/* Get the current model intensity */
		model_int = get_intensity(model_version);
167

168
169
170
171
		if ( pass == 1 ) {

			/* User asked for max only? */
			if ( !max_only ) {
172
173
				set_intensity(model_version,
				              model_int + intensity);
174
175
			} else {
				if ( intensity>get_intensity(model_version) ) {
176
					set_intensity(model_version, intensity);
177
				}
178
179
			}

180
181
182
183
			/* Increase redundancy */
			int cur_redundancy = get_redundancy(model_version);
			set_redundancy(model_version, cur_redundancy+1);

184

185
186
		} else if ( pass == 2 ) {

187
			double dev = get_temp1(model_version);
188
189
190
191
192
193
194
195
196

			/* Other ways of estimating the sigma are possible,
			 * choose from:
			 *    dev += pow(get_esd_intensity(refl), 2.0);
			 *    dev += pow(intensity, 2.0);
			 * But alter the other part of the calculation below
			 * as well. */
			dev += pow(intensity - model_int, 2.0);

197
			set_temp1(model_version, dev);
198
199
200
201
202
203
204
205
206

			if ( hist_vals != NULL ) {
				int p = *hist_n;
				if ( (h==hist_h) && (k==hist_k)
				  && (l==hist_l) )
				{
					hist_vals[p] = intensity;
					*hist_n = p+1;
				}
Thomas White's avatar
Thomas White committed
207

208
			}
209

210
211
		}

212
	}
213
214


215
216
217
}


218
219
220
enum {
	SCALE_NONE,
	SCALE_CONSTINT,
221
222
	SCALE_INTPERBRAGG,
	SCALE_TWOPASS,
223
224
225
};


Thomas White's avatar
Thomas White committed
226
static void scale_intensities(RefList *model, RefList *new, const SymOpList *sym)
227
228
{
	double s;
229
230
	double top = 0.0;
	double bot = 0.0;
231
	const int scaling = SCALE_INTPERBRAGG;
232
233
234
	Reflection *refl;
	RefListIterator *iter;
	Reflection *model_version;
235

236
237
238
	for ( refl = first_refl(new, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
239
240
241

		double i1, i2;
		signed int hu, ku, lu;
242
		signed int h, k, l;
243

244
		get_indices(refl, &h, &k, &l);
245

246
		switch ( scaling ) {
Thomas White's avatar
Thomas White committed
247
248

			case SCALE_TWOPASS :
249

250
251
252
			model_version = find_refl(model, h, k, l);
			if ( model_version == NULL ) continue;

Thomas White's avatar
Thomas White committed
253
			get_asymm(sym, h, k, l, &hu, &ku, &lu);
254

255
256
			i1 = get_intensity(model_version);
			i2 = get_intensity(refl);
257
258
259
260
261

			/* Calculate LSQ estimate of scaling factor */
			top += i1 * i2;
			bot += i2 * i2;
			break;
262

Thomas White's avatar
Thomas White committed
263
			case SCALE_CONSTINT :
264
			/* Sum up the intensity in the pattern */
265
			i2 = get_intensity(refl);
266
267
268
			top += i2;
			break;

Thomas White's avatar
Thomas White committed
269
			case SCALE_INTPERBRAGG :
270
			/* Sum up the intensity in the pattern */
271
			i2 = get_intensity(refl);
272
273
274
275
			top += i2;
			bot += 1.0;
			break;

276
		}
277
278
279

	}

280
	switch ( scaling ) {
Thomas White's avatar
Thomas White committed
281
282

		case SCALE_TWOPASS :
283
284
		s = top / bot;
		break;
Thomas White's avatar
Thomas White committed
285
286

		case SCALE_CONSTINT :
287
288
		s = 1000.0 / top;
		break;
Thomas White's avatar
Thomas White committed
289
290

		case SCALE_INTPERBRAGG :
291
292
		s = 1000.0 / (top/bot);
		break;
Thomas White's avatar
Thomas White committed
293

294
	}
295
296

	/* Multiply the new pattern up by "s" */
297
298
299
300
301
	for ( refl = first_refl(new, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {

		double intensity = get_intensity(refl);
302
		set_intensity(refl, intensity*s);
303

304
305
306
307
	}
}


308
static void merge_all(FILE *fh, RefList *model,
309
                      int config_maxonly, int config_scale, int config_sum,
310
                      int config_startafter, int config_stopafter,
Thomas White's avatar
Thomas White committed
311
                      const SymOpList *sym,
312
                      int n_total_patterns,
313
314
                      double *hist_vals, signed int hist_h,
                      signed int hist_k, signed int hist_l,
315
                      int *hist_i, int pass)
316
{
317
	int rval;
Thomas White's avatar
Thomas White committed
318
	int n_patterns = 0;
319
320
321
	int n_used = 0;
	Reflection *refl;
	RefListIterator *iter;
322

323
324
325
	if ( skip_some_files(fh, config_startafter) ) {
		ERROR("Failed to skip first %i files.\n", config_startafter);
		return;
326
327
	}

328
329
	do {

330
		struct image image;
331

332
333
		image.det = NULL;

334
335
		/* Get data from next chunk */
		rval = read_chunk(fh, &image);
Thomas White's avatar
Thomas White committed
336
		if ( rval ) break;
337

338
339
340
		n_patterns++;

		if ( (image.reflections != NULL) && (image.indexed_cell) ) {
341
342
343

			/* Scale if requested */
			if ( config_scale ) {
344
345
				scale_intensities(model, image.reflections,
				                  sym);
346
347
			}

348
349
			merge_pattern(model, image.reflections, config_maxonly,
			              sym, hist_vals, hist_h, hist_k, hist_l,
350
			              hist_i, pass);
351

352
			n_used++;
353
354
355

		}

Thomas White's avatar
Thomas White committed
356
		free(image.filename);
357
358
359
		reflist_free(image.reflections);
		image_feature_list_free(image.features);
		cell_free(image.indexed_cell);
360

361
362
		progress_bar(n_patterns, n_total_patterns-config_startafter,
		             "Merging");
363

364
	} while ( rval == 0 );
365

366
	/* Divide by counts to get mean intensity if necessary */
367
	if ( (pass == 1) && !config_sum && !config_maxonly ) {
368

369
370
		Reflection *refl;
		RefListIterator *iter;
371

372
373
374
		for ( refl = first_refl(model, &iter);
		      refl != NULL;
		      refl = next_refl(refl, iter) ) {
375

376
377
			double intensity = get_intensity(refl);
			int red = get_redundancy(refl);
378

379
			set_intensity(refl, intensity / (double)red);
380

Thomas White's avatar
Thomas White committed
381
		}
382

Thomas White's avatar
Thomas White committed
383
384
	}

385
	/* Calculate ESDs */
386
387
388
389
	if ( pass == 2 ) {
		for ( refl = first_refl(model, &iter);
		      refl != NULL;
		      refl = next_refl(refl, iter) ) {
390

391
			double sum_squared_dev = get_temp1(refl);
392
			int red = get_redundancy(refl);
393
			int h, k, l;
394
			double esd;
395
			get_indices(refl,&h,&k,&l);
396
397
398
399
400
401
402
403
404
405
406
407
408
409

			/* Other ways of estimating the sigma are possible,
			 * such as:
			 *
			 *    double intensity = get_intensity(refl);
			 *    esd = sqrt( (sum_squared_dev/(double)red)
			 *              - pow(intensity,2.0) ) );
			 *
			 * But alter the other part of the calculation above
			 * as well. */
			esd = sqrt(sum_squared_dev)/(double)red;

			set_esd_intensity(refl, esd);

410
		}
411
	}
Thomas White's avatar
Thomas White committed
412

413
414
415
	if ( pass == 1 ) {
		STATUS("%i of the patterns could be used.\n", n_used);
	}
416
417
418
}


419
420
421
422
int main(int argc, char *argv[])
{
	int c;
	char *filename = NULL;
423
	char *output = NULL;
424
	FILE *fh;
425
	RefList *model;
Thomas White's avatar
Thomas White committed
426
	int config_maxonly = 0;
427
	int config_startafter = 0;
Thomas White's avatar
Thomas White committed
428
	int config_stopafter = 0;
429
	int config_sum = 0;
430
	int config_scale = 0;
431
	unsigned int n_total_patterns;
Thomas White's avatar
Thomas White committed
432
433
	char *sym_str = NULL;
	SymOpList *sym;
434
	char *pdb = NULL;
435
436
	char *histo = NULL;
	signed int hist_h, hist_k, hist_l;
437
438
	signed int hist_nbins=50;
	float hist_min=0.0, hist_max=0.0;
439
440
	double *hist_vals = NULL;
	int hist_i;
441
	int space_for_hist = 0;
442
	char *histo_params = NULL;
443
444
445
446
447

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
448
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
449
		{"max-only",           0, &config_maxonly,     1},
Thomas White's avatar
Thomas White committed
450
		{"output-every",       1, NULL,               'e'},
Thomas White's avatar
Thomas White committed
451
		{"stop-after",         1, NULL,               's'},
452
		{"start-after",        1, NULL,               'f'},
453
		{"sum",                0, &config_sum,         1},
454
		{"scale",              0, &config_scale,       1},
455
		{"symmetry",           1, NULL,               'y'},
456
		{"histogram",          1, NULL,               'g'},
457
		{"hist-parameters",    1, NULL,               'z'},
458
459
460
461
		{0, 0, NULL, 0}
	};

	/* Short options */
462
	while ((c = getopt_long(argc, argv, "hi:e:o:p:y:g:f:b:z:",
463
	                        longopts, NULL)) != -1) {
464
465

		switch (c) {
Thomas White's avatar
Thomas White committed
466
		case 'h' :
467
468
469
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
470
		case 'i' :
471
472
473
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
474
		case 'o' :
475
476
477
			output = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
478
		case 's' :
Thomas White's avatar
Thomas White committed
479
480
481
			config_stopafter = atoi(optarg);
			break;

482
483
484
485
		case 'f' :
			config_startafter = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
486
		case 'p' :
487
488
489
			pdb = strdup(optarg);
			break;

490
		case 'y' :
Thomas White's avatar
Thomas White committed
491
			sym_str = strdup(optarg);
492
493
			break;

494
495
496
497
		case 'g' :
			histo = strdup(optarg);
			break;

498
499
500
501
		case 'z' :
			histo_params = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
502
		case 0 :
503
504
			break;

Thomas White's avatar
Thomas White committed
505
		default :
506
507
508
509
510
511
512
513
514
515
			return 1;
		}

	}

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

516
517
518
519
	if ( output == NULL ) {
		output = strdup("processed.hkl");
	}

Thomas White's avatar
Thomas White committed
520
521
522
	if ( sym_str == NULL ) sym_str = strdup("1");
	sym = get_pointgroup(sym_str);
	free(sym_str);
523

524
525
526
527
528
529
530
531
532
533
534
	/* Open the data stream */
	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
Tidy-up    
Thomas White committed
535

536
537
	/* Count the number of patterns in the file */
	n_total_patterns = count_patterns(fh);
Thomas White's avatar
Thomas White committed
538
539
540
541
	if ( n_total_patterns == 0 ) {
		ERROR("No patterns to process.\n");
		return 1;
	}
542
543
	STATUS("There are %i patterns to process\n", n_total_patterns);
	rewind(fh);
544

545
546
	model = reflist_new();

547
	if ( histo != NULL ) {
Thomas White's avatar
Thomas White committed
548

549
		int r;
Thomas White's avatar
Thomas White committed
550

551
552
553
554
555
		r = sscanf(histo, "%i,%i,%i", &hist_h, &hist_k, &hist_l);
		if ( r != 3 ) {
			ERROR("Invalid indices for '--histogram'\n");
			return 1;
		}
Thomas White's avatar
Thomas White committed
556

Thomas White's avatar
Thomas White committed
557
		space_for_hist = n_total_patterns * num_equivs(sym, NULL);
558
559
560
		hist_vals = malloc(space_for_hist * sizeof(double));
		free(histo);
		STATUS("Histogramming %i %i %i -> ", hist_h, hist_k, hist_l);
Thomas White's avatar
Thomas White committed
561

562
		/* Put into the asymmetric cell for the target group */
Thomas White's avatar
Thomas White committed
563
564
		get_asymm(sym, hist_h, hist_k, hist_l,
		          &hist_h, &hist_k, &hist_l);
565
		STATUS("%i %i %i\n", hist_h, hist_k, hist_l);
Thomas White's avatar
Thomas White committed
566

567
568
	}

569
	if ( histo_params != NULL ) {
Thomas White's avatar
Thomas White committed
570

571
		int rr;
Thomas White's avatar
Thomas White committed
572
573
574

		rr = sscanf(histo_params, "%f,%f,%i", &hist_min, &hist_max,
		                                      &hist_nbins);
575
576
577
578
579
		if ( rr != 3 ) {
			ERROR("Invalid parameters for '--hist-parameters'\n");
			return 1;
		}
		free(histo_params);
Thomas White's avatar
Thomas White committed
580
581
582
		if ( hist_max <= hist_min ) {
			ERROR("Invalid range for '--hist-parameters'. "
			      "Make sure that 'max' is greater than 'min'.\n");
583
584
			return 1;
		}
Thomas White's avatar
Thomas White committed
585

586
587
	}

588
	hist_i = 0;
589
	merge_all(fh, model, config_maxonly, config_scale, config_sum,
590
	          config_startafter, config_stopafter,
591
                  sym, n_total_patterns,
592
                  NULL, 0, 0, 0, NULL, 1);
Thomas White's avatar
Thomas White committed
593
594
595
596
	if ( ferror(fh) ) {
		ERROR("Stream read error.\n");
		return 1;
	}
597
	rewind(fh);
598

599
600
	STATUS("Extra pass to calculate ESDs...\n");
	rewind(fh);
601
602
	merge_all(fh, model, config_maxonly, config_scale, 0,
	          config_startafter, config_stopafter, sym, n_total_patterns,
603
	          hist_vals, hist_h, hist_k, hist_l, &hist_i, 2);
Thomas White's avatar
Thomas White committed
604
605
606
607
	if ( ferror(fh) ) {
		ERROR("Stream read error.\n");
		return 1;
	}
608

609
610
611
612
613
614
615
616
617
618
619
	if ( space_for_hist && (hist_i >= space_for_hist) ) {
		ERROR("Histogram array was too small!\n");
	}

	if ( hist_vals != NULL ) {
		STATUS("%i %i %i was seen %i times.\n", hist_h, hist_k, hist_l,
		                                        hist_i);
		plot_histogram(hist_vals, hist_i, hist_min, hist_max,
		               hist_nbins);
	}

620
	write_reflist(output, model);
621
622

	fclose(fh);
623

624
	free(sym);
625
	reflist_free(model);
Thomas White's avatar
Thomas White committed
626
	free(output);
627

628
629
	return 0;
}