partialator.c 13.3 KB
Newer Older
1
/*
2
 * partialator.c
3
 *
4
 * Scaling and post refinement for coherent nanocrystallography
5
 *
6
7
8
9
 * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
 *                  a research centre of the Helmholtz Association.
 *
 * Authors:
10
 *   2010-2012 Thomas White <taw@physics.org>
11
 *
Thomas White's avatar
Thomas White committed
12
13
14
15
16
17
18
19
20
21
22
23
24
25
 * 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/>.
26
27
28
29
30
31
32
33
34
35
36
37
38
39
 *
 */


#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>
Thomas White's avatar
Thomas White committed
40
#include <assert.h>
Thomas White's avatar
Tidy up    
Thomas White committed
41
#include <pthread.h>
42
#include <gsl/gsl_errno.h>
43

Thomas White's avatar
Thomas White committed
44
45
46
47
48
49
50
51
52
53
54
#include <utils.h>
#include <hdf5-file.h>
#include <symmetry.h>
#include <stream.h>
#include <geometry.h>
#include <peaks.h>
#include <thread-pool.h>
#include <beam-parameters.h>
#include <reflist.h>
#include <reflist-utils.h>

55
#include "post-refinement.h"
56
#include "hrs-scaling.h"
Thomas White's avatar
Thomas White committed
57
#include "scaling-report.h"
58
59
60
61
62
63


static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
64
"Scaling and post refinement for coherent nanocrystallography.\n"
65
66
67
"\n"
"  -h, --help                 Display this help message.\n"
"\n"
Thomas White's avatar
Thomas White committed
68
69
"  -i, --input=<filename>     Specify the name of the input 'stream'.\n"
"                              (must be a file, not e.g. stdin)\n"
Thomas White's avatar
Thomas White committed
70
"  -o, --output=<filename>    Output filename.  Default: partialator.hkl.\n"
Thomas White's avatar
Thomas White committed
71
"  -g. --geometry=<file>      Get detector geometry from file.\n"
72
73
74
75
"  -b, --beam=<file>          Get beam parameters from file, which provides\n"
"                              initial values for parameters, and nominal\n"
"                              wavelengths if no per-shot value is found in \n"
"                              an HDF5 file.\n"
76
"  -y, --symmetry=<sym>       Merge according to symmetry <sym>.\n"
77
"  -n, --iterations=<n>       Run <n> cycles of scaling and post-refinement.\n"
78
"      --no-scale             Fix all the scaling factors at unity.\n"
79
"  -r, --reference=<file>     Refine images against reflections in <file>,\n"
80
81
"                              instead of taking the mean of the intensity\n"
"                              estimates.\n"
82
"\n"
Thomas White's avatar
Thomas White committed
83
84
85
86
"  -j <n>                     Run <n> analyses in parallel.\n");
}


Thomas White's avatar
Thomas White committed
87
struct refine_args
Thomas White's avatar
Thomas White committed
88
{
89
	RefList *full;
Thomas White's avatar
Thomas White committed
90
91
92
93
	struct image *image;
};


Thomas White's avatar
Thomas White committed
94
struct queue_args
Thomas White's avatar
Thomas White committed
95
{
Thomas White's avatar
Thomas White committed
96
97
98
99
100
101
102
103
104
105
106
	int n;
	int n_done;
	int n_total_patterns;
	struct image *images;
	struct refine_args task_defaults;
};


static void refine_image(void *task, int id)
{
	struct refine_args *pargs = task;
107
	struct image *image = pargs->image;
Thomas White's avatar
Thomas White committed
108
	image->id = id;
109

110
	pr_refine(image, pargs->full);
111
112
113
}


Thomas White's avatar
Thomas White committed
114
static void *get_image(void *vqargs)
115
{
Thomas White's avatar
Thomas White committed
116
117
	struct refine_args *task;
	struct queue_args *qargs = vqargs;
Thomas White's avatar
Thomas White committed
118

Thomas White's avatar
Thomas White committed
119
120
	task = malloc(sizeof(struct refine_args));
	memcpy(task, &qargs->task_defaults, sizeof(struct refine_args));
121

Thomas White's avatar
Thomas White committed
122
	task->image = &qargs->images[qargs->n];
123

Thomas White's avatar
Thomas White committed
124
125
126
127
	qargs->n++;

	return task;
}
Thomas White's avatar
Thomas White committed
128

129

Thomas White's avatar
Thomas White committed
130
131
132
133
134
135
static void done_image(void *vqargs, void *task)
{
	struct queue_args *qargs = vqargs;

	qargs->n_done++;

Thomas White's avatar
Thomas White committed
136
	progress_bar(qargs->n_done, qargs->n_total_patterns, "Refining");
Thomas White's avatar
Thomas White committed
137
138
139
140
141
	free(task);
}


static void refine_all(struct image *images, int n_total_patterns,
142
                       struct detector *det,
Thomas White's avatar
Thomas White committed
143
                       RefList *full, int nthreads)
Thomas White's avatar
Thomas White committed
144
145
146
147
148
149
150
151
152
153
{
	struct refine_args task_defaults;
	struct queue_args qargs;

	task_defaults.full = full;
	task_defaults.image = NULL;

	qargs.task_defaults = task_defaults;
	qargs.n = 0;
	qargs.n_done = 0;
Thomas White's avatar
Thomas White committed
154
155
	qargs.n_total_patterns = n_total_patterns;
	qargs.images = images;
Thomas White's avatar
Thomas White committed
156

Thomas White's avatar
Thomas White committed
157
	/* Don't have threads which are doing nothing */
Thomas White's avatar
Whoops    
Thomas White committed
158
	if ( n_total_patterns < nthreads ) nthreads = n_total_patterns;
Thomas White's avatar
Thomas White committed
159

Thomas White's avatar
Thomas White committed
160
	run_threads(nthreads, refine_image, get_image, done_image,
Thomas White's avatar
Thomas White committed
161
	            &qargs, n_total_patterns, 0, 0, 0);
162
163
164
}


165
/* Decide which reflections can be scaled */
166
static int select_scalable_reflections(RefList *list, RefList *reference)
167
168
169
170
171
172
173
174
175
{
	Reflection *refl;
	RefListIterator *iter;
	int nobs = 0;

	for ( refl = first_refl(list, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {

176
		int sc = 1;
177
		double v, esd;
178

179
180
181
		/* This means the reflection was not found on the last check */
		if ( get_redundancy(refl) == 0 ) sc = 0;

182
		if ( get_partiality(refl) < 0.1 ) sc = 0;
183
		v = fabs(get_intensity(refl));
184
185
		esd = get_esd_intensity(refl);
		if ( v < 0.5*esd ) sc = 0;
186

187
188
189
		/* If we are scaling against a reference set, we additionally
		 * require that this reflection is in the reference list. */
		if ( reference != NULL ) {
190
191
			signed int h, k, l;
			get_indices(refl, &h, &k, &l);
192
			if ( find_refl(reference, h, k, l) == NULL ) sc = 0;
193
194
		}

195
196
197
		set_scalable(refl, sc);

		if ( sc ) nobs++;
198
199
200
201
202
203
	}

	return nobs;
}


204
static void select_reflections_for_refinement(struct image *images, int n,
205
                                              RefList *full, int have_reference)
206
207
208
209
210
211
212
{
	int i;

	for ( i=0; i<n; i++ ) {

		Reflection *refl;
		RefListIterator *iter;
213
214
215
216
217
		int n_acc = 0;
		int n_nomatch = 0;
		int n_noscale = 0;
		int n_fewmatch = 0;
		int n_ref = 0;
218

219
220
		if ( images[i].pr_dud ) continue;

221
222
223
224
225
226
227
		for ( refl = first_refl(images[i].reflections, &iter);
		      refl != NULL;
		      refl = next_refl(refl, iter) )
		{
			signed int h, k, l;
			int sc;

228
229
230
231
232
233
			n_ref++;

			/* We require that the reflection itself is scalable
			 * (i.e. sensible partiality and intensity) and that
			 * the "full" estimate of this reflection is made from
			 * at least two parts. */
234
235
			get_indices(refl, &h, &k, &l);
			sc = get_scalable(refl);
236
237
238
			if ( !sc ) {

				n_noscale++;
239
				set_refinable(refl, 0);
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256

			} else {

				Reflection *f = find_refl(full, h, k, l);

				if ( f != NULL ) {

					int r = get_redundancy(f);
					if ( (r >= 2) || have_reference ) {
						set_refinable(refl, 1);
						n_acc++;
					} else {
						n_fewmatch++;
					}

				} else {
					n_nomatch++;
257
					set_refinable(refl, 0);
258
				}
259
260
261
262

			}
		}

263
264
265
		//STATUS("Image %4i: %i guide reflections accepted "
		//       "(%i not scalable, %i few matches, %i total)\n",
		//       i, n_acc, n_noscale, n_fewmatch, n_ref);
266
267
268
269
270

		/* This would be a silly situation, since there must be a match
		 * if THIS pattern has a scalable part of the reflection! */
		assert(n_nomatch == 0);

271
272
273
274
	}
}


Thomas White's avatar
Thomas White committed
275
276
277
278
int main(int argc, char *argv[])
{
	int c;
	char *infile = NULL;
Thomas White's avatar
Thomas White committed
279
	char *outfile = NULL;
Thomas White's avatar
Thomas White committed
280
	char *geomfile = NULL;
Thomas White's avatar
Thomas White committed
281
282
	char *sym_str = NULL;
	SymOpList *sym;
Thomas White's avatar
Thomas White committed
283
	FILE *fh;
Thomas White's avatar
Thomas White committed
284
285
	int nthreads = 1;
	struct detector *det;
Thomas White's avatar
Thomas White committed
286
287
	int i;
	int n_total_patterns;
288
289
	struct image *images;
	int n_iter = 10;
Thomas White's avatar
Thomas White committed
290
	struct beam_params *beam = NULL;
291
	RefList *full;
Thomas White's avatar
Thomas White committed
292
	int n_usable_patterns = 0;
293
	int nobs;
294
	char *reference_file = NULL;
295
	RefList *reference = NULL;
296
	int n_dud;
297
	int have_reference = 0;
Thomas White's avatar
Thomas White committed
298
299
	char cmdline[1024];
	SRContext *sr;
300
	int noscale = 0;
Thomas White's avatar
Thomas White committed
301
302
303
304
305

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
306
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
307
		{"geometry",           1, NULL,               'g'},
Thomas White's avatar
Thomas White committed
308
		{"beam",               1, NULL,               'b'},
309
		{"symmetry",           1, NULL,               'y'},
310
		{"iterations",         1, NULL,               'n'},
311
		{"no-scale",           0, &noscale,            1},
312
		{"reference",          1, NULL,               'r'},
Thomas White's avatar
Thomas White committed
313
314
315
		{0, 0, NULL, 0}
	};

Thomas White's avatar
Thomas White committed
316
317
318
319
320
321
	cmdline[0] = '\0';
	for ( i=1; i<argc; i++ ) {
		strncat(cmdline, argv[i], 1023-strlen(cmdline));
		strncat(cmdline, " ", 1023-strlen(cmdline));
	}

Thomas White's avatar
Thomas White committed
322
	/* Short options */
Thomas White's avatar
Thomas White committed
323
	while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:r:j:",
324
325
	                        longopts, NULL)) != -1)
	{
Thomas White's avatar
Thomas White committed
326
327

		switch (c) {
Thomas White's avatar
Thomas White committed
328
329

			case 'h' :
Thomas White's avatar
Thomas White committed
330
331
332
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
333
			case 'i' :
Thomas White's avatar
Thomas White committed
334
335
336
			infile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
337
			case 'g' :
Thomas White's avatar
Thomas White committed
338
339
340
			geomfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
341
			case 'j' :
Thomas White's avatar
Thomas White committed
342
343
344
			nthreads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
345
			case 'y' :
Thomas White's avatar
Thomas White committed
346
			sym_str = strdup(optarg);
347
348
			break;

Thomas White's avatar
Thomas White committed
349
			case 'o' :
Thomas White's avatar
Thomas White committed
350
351
352
			outfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
353
			case 'n' :
354
355
356
			n_iter = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
357
			case 'b' :
Thomas White's avatar
Thomas White committed
358
359
360
361
362
363
364
365
			beam = get_beam_parameters(optarg);
			if ( beam == NULL ) {
				ERROR("Failed to load beam parameters"
				      " from '%s'\n", optarg);
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
366
			case 'r' :
367
			reference_file = strdup(optarg);
368
369
			break;

Thomas White's avatar
Thomas White committed
370
			case 0 :
Thomas White's avatar
Thomas White committed
371
372
			break;

373
374
375
			case '?' :
			break;

Thomas White's avatar
Thomas White committed
376
			default :
377
378
			ERROR("Unhandled option '%c'\n", c);
			break;
Thomas White's avatar
Thomas White committed
379

Thomas White's avatar
Thomas White committed
380
381
382
383
		}

	}

Thomas White's avatar
Thomas White committed
384
385
386
387
388
	if ( nthreads < 1 ) {
		ERROR("Invalid number of threads.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
389
	/* Sanitise input filename and open */
Thomas White's avatar
Thomas White committed
390
391
392
393
394
395
396
397
398
399
400
401
402
	if ( infile == NULL ) {
		infile = strdup("-");
	}
	if ( strcmp(infile, "-") == 0 ) {
		fh = stdin;
	} else {
		fh = fopen(infile, "r");
	}
	if ( fh == NULL ) {
		ERROR("Failed to open input file '%s'\n", infile);
		return 1;
	}

Thomas White's avatar
Thomas White committed
403
404
	/* Sanitise output filename */
	if ( outfile == NULL ) {
Thomas White's avatar
Whoops    
Thomas White committed
405
		outfile = strdup("partialator.hkl");
Thomas White's avatar
Thomas White committed
406
407
	}

Thomas White's avatar
Thomas White committed
408
409
410
	if ( sym_str == NULL ) sym_str = strdup("1");
	sym = get_pointgroup(sym_str);
	free(sym_str);
Thomas White's avatar
Thomas White committed
411

Thomas White's avatar
Thomas White committed
412
	/* Get detector geometry */
Thomas White's avatar
Thomas White committed
413
414
415
416
417
418
419
	det = get_detector_geometry(geomfile);
	if ( det == NULL ) {
		ERROR("Failed to read detector geometry from '%s'\n", geomfile);
		return 1;
	}
	free(geomfile);

Thomas White's avatar
Thomas White committed
420
421
422
423
424
	if ( beam == NULL ) {
		ERROR("You must provide a beam parameters file.\n");
		return 1;
	}

425
	if ( reference_file != NULL ) {
Thomas White's avatar
Thomas White committed
426

427
		RefList *list;
428

429
		list = read_reflections(reference_file);
430
431
432
433
		if ( list == NULL ) {
			ERROR("Failed to read '%s'\n", reference_file);
			return 1;
		}
434
		free(reference_file);
435
		reference = asymmetric_indices(list, sym);
436
		reflist_free(list);
437
		have_reference = 1;
438

439
440
	}

Thomas White's avatar
Thomas White committed
441
	n_total_patterns = count_patterns(fh);
Thomas White's avatar
Thomas White committed
442
443
444
445
	if ( n_total_patterns == 0 ) {
		ERROR("No patterns to process.\n");
		return 1;
	}
Thomas White's avatar
Thomas White committed
446
447
	STATUS("There are %i patterns to process\n", n_total_patterns);

448
449
	gsl_set_error_handler_off();

450
451
452
453
454
455
456
457
	images = malloc(n_total_patterns * sizeof(struct image));
	if ( images == NULL ) {
		ERROR("Couldn't allocate memory for images.\n");
		return 1;
	}

	/* Fill in what we know about the images so far */
	rewind(fh);
458
	nobs = 0;
459
460
	for ( i=0; i<n_total_patterns; i++ ) {

461
		RefList *as;
Thomas White's avatar
Thomas White committed
462
		struct image *cur;
463

Thomas White's avatar
Thomas White committed
464
		cur = &images[n_usable_patterns];
465

Thomas White's avatar
Thomas White committed
466
467
468
		cur->det = det;

		if ( read_chunk(fh, cur) != 0 ) {
Thomas White's avatar
Thomas White committed
469
470
			/* Should not happen, because we counted the patterns
			 * earlier. */
Thomas White's avatar
Thomas White committed
471
			ERROR("Failed to read chunk from the input stream.\n");
472
473
474
			return 1;
		}

475
		/* Won't be needing this, if it exists */
Thomas White's avatar
Thomas White committed
476
477
		image_feature_list_free(cur->features);
		cur->features = NULL;
478

Thomas White's avatar
Thomas White committed
479
		/* "n_usable_patterns" will not be incremented in this case */
Thomas White's avatar
Thomas White committed
480
		if ( cur->indexed_cell == NULL ) continue;
481

Thomas White's avatar
Thomas White committed
482
		/* Fill in initial estimates of stuff */
Thomas White's avatar
Thomas White committed
483
484
485
486
487
		cur->div = beam->divergence;
		cur->bw = beam->bandwidth;
		cur->width = det->max_fs;
		cur->height = det->max_ss;
		cur->osf = 1.0;
488
		cur->profile_radius = beam->profile_radius;
Thomas White's avatar
Thomas White committed
489
		cur->pr_dud = 0;
490

Thomas White's avatar
Thomas White committed
491
		/* Muppet proofing */
Thomas White's avatar
Thomas White committed
492
493
494
		cur->data = NULL;
		cur->flags = NULL;
		cur->beam = NULL;
495

Thomas White's avatar
Thomas White committed
496
		/* This is the raw list of reflections */
Thomas White's avatar
Thomas White committed
497
498
499
		as = asymmetric_indices(cur->reflections, sym);
		reflist_free(cur->reflections);
		cur->reflections = as;
500

501
		update_partialities(cur);
502

503
504
		nobs += select_scalable_reflections(cur->reflections,
		                                    reference);
505

506
		progress_bar(i, n_total_patterns-1, "Loading pattern data");
Thomas White's avatar
Thomas White committed
507
		n_usable_patterns++;
508
509
510

	}
	fclose(fh);
Thomas White's avatar
Thomas White committed
511

512
	/* Make initial estimates */
513
	STATUS("Performing initial scaling.\n");
514
	if ( noscale ) STATUS("Scale factors fixed at 1.\n");
Thomas White's avatar
Thomas White committed
515
	full = scale_intensities(images, n_usable_patterns, reference,
516
	                         nthreads, noscale);
Thomas White's avatar
Thomas White committed
517

518
519
520
	sr = sr_titlepage(images, n_usable_patterns, "scaling-report.pdf",
	                  infile, cmdline);
	sr_iteration(sr, 0, images, n_usable_patterns, full);
Thomas White's avatar
Thomas White committed
521

Thomas White's avatar
Thomas White committed
522
	/* Iterate */
523
	for ( i=0; i<n_iter; i++ ) {
Thomas White's avatar
Thomas White committed
524

525
		int j;
526
		RefList *comp;
527

528
		STATUS("Post refinement cycle %i of %i\n", i+1, n_iter);
Thomas White's avatar
Thomas White committed
529

530
531
532
533
534
		if ( reference == NULL ) {
			comp = full;
		} else {
			comp = reference;
		}
535

Thomas White's avatar
Thomas White committed
536
		/* Refine the geometry of all patterns to get the best fit */
537
538
		select_reflections_for_refinement(images, n_usable_patterns,
		                                  comp, have_reference);
Thomas White's avatar
Thomas White committed
539
		refine_all(images, n_usable_patterns, det, comp, nthreads);
Thomas White's avatar
Thomas White committed
540

541
		nobs = 0;
542
543
544
545
		for ( j=0; j<n_usable_patterns; j++ ) {

			struct image *cur = &images[j];

546
			nobs += select_scalable_reflections(cur->reflections,
547
			                                    reference);
548
549
550

		}

Thomas White's avatar
Thomas White committed
551
		/* Re-estimate all the full intensities */
552
		reflist_free(full);
Thomas White's avatar
Thomas White committed
553
		full = scale_intensities(images, n_usable_patterns,
554
		                         reference, nthreads, noscale);
555

556
557
		sr_iteration(sr, i+1, images, n_usable_patterns, full);

Thomas White's avatar
Thomas White committed
558
559
	}

560
561
	sr_finish(sr);

562
	n_dud = 0;
Thomas White's avatar
Thomas White committed
563
	for ( i=0; i<n_usable_patterns; i++ ) {
564
		if ( images[i].pr_dud ) n_dud++;
Thomas White's avatar
Thomas White committed
565
	}
566
	STATUS("%i images could not be refined on the last cycle.\n", n_dud);
Thomas White's avatar
Thomas White committed
567
568

	/* Output results */
569
	write_reflist(outfile, full);
Thomas White's avatar
Thomas White committed
570

Thomas White's avatar
Thomas White committed
571
	/* Clean up */
Thomas White's avatar
Thomas White committed
572
	for ( i=0; i<n_usable_patterns; i++ ) {
Thomas White's avatar
Thomas White committed
573
		reflist_free(images[i].reflections);
Thomas White's avatar
Thomas White committed
574
	}
575
	reflist_free(full);
Thomas White's avatar
Thomas White committed
576
	free(sym);
Thomas White's avatar
Thomas White committed
577
	free(outfile);
Thomas White's avatar
Thomas White committed
578
	free_detector_geometry(det);
579
	free(beam);
580
	if ( reference != NULL ) {
581
		reflist_free(reference);
582
	}
Thomas White's avatar
Thomas White committed
583
	for ( i=0; i<n_usable_patterns; i++ ) {
584
585
586
587
		cell_free(images[i].indexed_cell);
		free(images[i].filename);
	}
	free(images);
Thomas White's avatar
Thomas White committed
588
	free(infile);
589
590
591

	return 0;
}