partial_sim.c 13.4 KB
Newer Older
1
2
3
4
5
/*
 * partial_sim.c
 *
 * Generate partials for testing scaling
 *
6
7
8
9
 * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
 *                  a research centre of the Helmholtz Association.
 *
 * Authors:
10
 *   2011-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
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 <assert.h>
41
#include <pthread.h>
42

Thomas White's avatar
Thomas White committed
43
44
45
46
47
48
49
50
#include <utils.h>
#include <reflist-utils.h>
#include <symmetry.h>
#include <beam-parameters.h>
#include <detector.h>
#include <geometry.h>
#include <stream.h>
#include <thread-pool.h>
51

52
53
54
/* Number of bins for partiality graph */
#define NBINS 50

55

56
static void mess_up_cell(UnitCell *cell, double cnoise)
57
58
59
60
61
{
	double ax, ay, az;
	double bx, by, bz;
	double cx, cy, cz;

Thomas White's avatar
Thomas White committed
62
63
	//STATUS("Real:\n");
	//cell_print(cell);
Thomas White's avatar
Thomas White committed
64

65
	cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
66
67
68
69
70
71
72
73
74
	ax = flat_noise(ax, cnoise*fabs(ax)/100.0);
	ay = flat_noise(ay, cnoise*fabs(ay)/100.0);
	az = flat_noise(az, cnoise*fabs(az)/100.0);
	bx = flat_noise(bx, cnoise*fabs(bx)/100.0);
	by = flat_noise(by, cnoise*fabs(by)/100.0);
	bz = flat_noise(bz, cnoise*fabs(bz)/100.0);
	cx = flat_noise(cx, cnoise*fabs(cx)/100.0);
	cy = flat_noise(cy, cnoise*fabs(cy)/100.0);
	cz = flat_noise(cz, cnoise*fabs(cz)/100.0);
75
	cell_set_reciprocal(cell, ax, ay, az, bx, by, bz, cx, cy, cz);
Thomas White's avatar
Thomas White committed
76

Thomas White's avatar
Thomas White committed
77
78
	//STATUS("Changed:\n");
	//cell_print(cell);
79
80
}

81

82
83
/* For each reflection in "partial", fill in what the intensity would be
 * according to "full" */
84
static void calculate_partials(RefList *partial, double osf,
Thomas White's avatar
Thomas White committed
85
                               RefList *full, const SymOpList *sym,
86
                               int random_intensities,
87
88
                               pthread_mutex_t *full_lock,
                               unsigned long int *n_ref, double *p_hist,
89
                               double *p_max, double max_q, UnitCell *cell)
90
91
92
93
94
95
{
	Reflection *refl;
	RefListIterator *iter;

	for ( refl = first_refl(partial, &iter);
	      refl != NULL;
Thomas White's avatar
Thomas White committed
96
97
	      refl = next_refl(refl, iter) )
	{
98
99
		signed int h, k, l;
		Reflection *rfull;
Thomas White's avatar
Thomas White committed
100
		double p, Ip, If;
101
		int bin;
102
103

		get_indices(refl, &h, &k, &l);
Thomas White's avatar
Thomas White committed
104
		get_asymm(sym, h, k, l, &h, &k, &l);
105
106
		p = get_partiality(refl);

107
		pthread_mutex_lock(full_lock);
108
		rfull = find_refl(full, h, k, l);
109
110
		pthread_mutex_unlock(full_lock);

111
		if ( rfull == NULL ) {
112
			if ( random_intensities ) {
113
114
115
116
117
118

				/* The full reflection is immutable (in this
				 * program) once created, but creating it must
				 * be an atomic operation.  So do the whole
				 * thing under lock. */
				pthread_mutex_lock(full_lock);
119
120
				rfull = add_refl(full, h, k, l);
				If = fabs(gaussian_noise(0.0, 1000.0));
121
				set_intensity(rfull, If);
122
				set_redundancy(rfull, 1);
123
124
				pthread_mutex_unlock(full_lock);

125
126
127
128
			} else {
				set_redundancy(refl, 0);
				If = 0.0;
			}
129
130
		} else {
			If = get_intensity(rfull);
131
			if ( random_intensities ) {
132
				lock_reflection(rfull);
133
134
				int red = get_redundancy(rfull);
				set_redundancy(rfull, red+1);
135
				unlock_reflection(rfull);
136
			}
137
		}
138
139

		Ip = osf * p * If;
140

141
142
143
144
		bin = NBINS*2.0*resolution(cell, h, k, l) / max_q;
		if ( (bin < NBINS) && (bin>=0) ) {
			p_hist[bin] += p;
			n_ref[bin]++;
145
			if ( p > p_max[bin] ) p_max[bin] = p;
Thomas White's avatar
Thomas White committed
146
147
148
		} else {
			STATUS("Reflection out of histogram range: %e %i %f\n",
			       resolution(cell, h, k, l), bin,  p);
149
150
		}

151
		Ip = gaussian_noise(Ip, 100.0);
152

153
		set_intensity(refl, Ip);
154
		set_esd_intensity(refl, 100.0);
155
156
157
158
159
160
161
162
163
164
165
166
167
	}
}


static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
"Generate a stream containing partials from a reflection list.\n"
"\n"
" -h, --help              Display this help message.\n"
"\n"
"You need to provide the following basic options:\n"
168
169
170
171
172
173
" -i, --input=<file>       Read reflections from <file>.\n"
"                           Default: generate random ones instead (see -r).\n"
" -o, --output=<file>      Write partials in stream format to <file>.\n"
" -g. --geometry=<file>    Get detector geometry from file.\n"
" -b, --beam=<file>        Get beam parameters from file\n"
" -p, --pdb=<file>         PDB file from which to get the unit cell.\n"
174
"\n"
175
176
177
" -y, --symmetry=<sym>     Symmetry of the input reflection list.\n"
" -n <n>                   Simulate <n> patterns.  Default: 2.\n"
" -r, --save-random=<file> Save randomly generated intensities to file.\n"
178
179
180
" -c, --cnoise=<val>       Add random noise, with a flat distribution, to the\n"
"                          reciprocal lattice vector components given in the\n"
"                          stream, with maximum error +/- <val> percent.\n"
Thomas White's avatar
Thomas White committed
181
182
"     --pgraph=<file>      Write reflection counts and partiality statistics\n"
"                           to <file>.\n"
183
"\n"
184
185
186
187
);
}


188
189
190
191
192
193
struct queue_args
{
	RefList *full;
	pthread_mutex_t full_lock;

	int n_done;
194
	int n_started;
195
196
197
198
199
	int n_to_do;

	SymOpList *sym;
	int random_intensities;
	UnitCell *cell;
200
	double cnoise;
201
202

	struct image *template_image;
203
204
205
206
207
	double max_q;

	/* The overall histogram */
	double p_hist[NBINS];
	unsigned long int n_ref[NBINS];
208
	double p_max[NBINS];
209
210
211
212
213
214
215
216
217

	FILE *stream;
};


struct worker_args
{
	struct queue_args *qargs;
	struct image image;
218
219
220
221

	/* Histogram for this image */
	double p_hist[NBINS];
	unsigned long int n_ref[NBINS];
222
	double p_max[NBINS];
223
224
225
226
227
228
229
230
};


static void *create_job(void *vqargs)
{
	struct worker_args *wargs;
	struct queue_args *qargs = vqargs;

231
232
233
	/* All done already? */
	if ( qargs->n_started == qargs->n_to_do ) return NULL;

234
235
236
237
238
	wargs = malloc(sizeof(struct worker_args));

	wargs->qargs = qargs;
	wargs->image = *qargs->template_image;

239
240
	qargs->n_started++;

241
242
243
244
245
246
247
248
249
250
	return wargs;
}


static void run_job(void *vwargs, int cookie)
{
	double osf;
	struct quaternion orientation;
	struct worker_args *wargs = vwargs;
	struct queue_args *qargs = wargs->qargs;
251
	int i;
252
253
254
255
256
257
258
259
260
261

	osf = gaussian_noise(1.0, 0.3);

	/* Set up a random orientation */
	orientation = random_quaternion();
	wargs->image.indexed_cell = cell_rotate(qargs->cell, orientation);

	snprintf(wargs->image.filename, 255, "dummy.h5");
	wargs->image.reflections = find_intersections(&wargs->image,
	                                       wargs->image.indexed_cell);
262
263
264
265

	for ( i=0; i<NBINS; i++ ) {
		wargs->n_ref[i] = 0;
		wargs->p_hist[i] = 0.0;
266
		wargs->p_max[i] = 0.0;
267
268
	}

269
270
	calculate_partials(wargs->image.reflections, osf, qargs->full,
	                   qargs->sym, qargs->random_intensities,
271
	                   &qargs->full_lock,
272
273
	                   wargs->n_ref, wargs->p_hist, wargs->p_max,
	                   qargs->max_q, wargs->image.indexed_cell);
274
275

	/* Give a slightly incorrect cell in the stream */
276
	mess_up_cell(wargs->image.indexed_cell, qargs->cnoise);
277
278
279
280
281
282
283
}


static void finalise_job(void *vqargs, void *vwargs)
{
	struct worker_args *wargs = vwargs;
	struct queue_args *qargs = vqargs;
284
	int i;
285

286
	write_chunk(qargs->stream, &wargs->image, NULL, STREAM_INTEGRATED);
287

288
289
290
	for ( i=0; i<NBINS; i++ ) {
		qargs->n_ref[i] += wargs->n_ref[i];
		qargs->p_hist[i] += wargs->p_hist[i];
291
292
293
		if ( wargs->p_max[i] > qargs->p_max[i] ) {
			qargs->p_max[i] = wargs->p_max[i];
		}
294
295
	}

296
297
	qargs->n_done++;
	progress_bar(qargs->n_done, qargs->n_to_do, "Simulating");
298
299
300
301

	reflist_free(wargs->image.reflections);
	cell_free(wargs->image.indexed_cell);
	free(wargs);
302
303
304
}


305
306
307
308
309
310
311
312
313
314
int main(int argc, char *argv[])
{
	int c;
	char *input_file = NULL;
	char *output_file = NULL;
	char *beamfile = NULL;
	char *geomfile = NULL;
	char *cellfile = NULL;
	struct detector *det = NULL;
	struct beam_params *beam = NULL;
315
	RefList *full = NULL;
Thomas White's avatar
Thomas White committed
316
317
	char *sym_str = NULL;
	SymOpList *sym;
318
319
	UnitCell *cell = NULL;
	FILE *ofh;
320
	int n = 2;
321
	int random_intensities = 0;
322
	char *save_file = NULL;
323
324
325
	struct queue_args qargs;
	struct image image;
	int n_threads = 1;
326
327
	double cnoise = 0.0;
	char *rval;
328
329
330
	int i;
	FILE *fh;
	char *phist_file = NULL;
331
332
333
334
335
336
337
338
339
340

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"output",             1, NULL,               'o'},
		{"input",              1, NULL,               'i'},
		{"beam",               1, NULL,               'b'},
		{"pdb",                1, NULL,               'p'},
		{"geometry",           1, NULL,               'g'},
		{"symmetry",           1, NULL,               'y'},
341
		{"save-random",        1, NULL,               'r'},
342
		{"pgraph",             1, NULL,                2},
343
		{"cnoise",             1, NULL,               'c'},
344
345
346
347
		{0, 0, NULL, 0}
	};

	/* Short options */
348
	while ((c = getopt_long(argc, argv, "hi:o:b:p:g:y:n:r:j:c:",
Thomas White's avatar
Thomas White committed
349
350
	                        longopts, NULL)) != -1)
	{
351
		switch (c) {
Thomas White's avatar
Thomas White committed
352
353

			case 'h' :
354
355
356
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
357
			case 'o' :
358
359
360
			output_file = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
361
			case 'i' :
362
363
364
			input_file = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
365
			case 'b' :
366
367
368
			beamfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
369
			case 'p' :
370
371
372
			cellfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
373
			case 'g' :
374
375
376
			geomfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
377
			case 'y' :
Thomas White's avatar
Thomas White committed
378
			sym_str = strdup(optarg);
379
380
			break;

Thomas White's avatar
Thomas White committed
381
			case 'n' :
382
383
384
			n = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
385
			case 'r' :
386
387
388
			save_file = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
389
			case 'j' :
390
391
392
			n_threads = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
393
			case 'c' :
394
395
396
397
398
399
400
			cnoise = strtod(optarg, &rval);
			if ( *rval != '\0' ) {
				ERROR("Invalid cell noise value.\n");
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
401
			case 2 :
402
403
404
			phist_file = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
405
			case 0 :
406
407
			break;

Thomas White's avatar
Thomas White committed
408
			default :
409
			return 1;
Thomas White's avatar
Thomas White committed
410

411
412
413
		}
	}

414
415
416
417
418
	if ( n_threads < 1 ) {
		ERROR("Invalid number of threads.\n");
		return 1;
	}

419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
	/* Load beam */
	if ( beamfile == NULL ) {
		ERROR("You need to provide a beam parameters file.\n");
		return 1;
	}
	beam = get_beam_parameters(beamfile);
	if ( beam == NULL ) {
		ERROR("Failed to load beam parameters from '%s'\n", beamfile);
		return 1;
	}
	free(beamfile);

	/* Load cell */
	if ( cellfile == NULL ) {
		ERROR("You need to give a PDB file with the unit cell.\n");
		return 1;
	}
	cell = load_cell_from_pdb(cellfile);
	if ( cell == NULL ) {
		ERROR("Failed to get cell from '%s'\n", cellfile);
		return 1;
	}
	free(cellfile);

Thomas White's avatar
Thomas White committed
443
444
445
446
447
448
	if ( !cell_is_sensible(cell) ) {
		ERROR("Invalid unit cell parameters:\n");
		cell_print(cell);
		return 1;
	}

449
450
451
452
453
454
	/* Load geometry */
	if ( geomfile == NULL ) {
		ERROR("You need to give a geometry file.\n");
		return 1;
	}
	det = get_detector_geometry(geomfile);
455
	if ( det == NULL ) {
456
457
458
459
460
		ERROR("Failed to read geometry from '%s'\n", geomfile);
		return 1;
	}
	free(geomfile);

Thomas White's avatar
Thomas White committed
461
462
463
464
	if ( sym_str == NULL ) sym_str = strdup("1");
	sym = get_pointgroup(sym_str);
	free(sym_str);

465
466
	if ( save_file == NULL ) save_file = strdup("partial_sim.hkl");

467
	/* Load (full) reflections */
468
469
470
471
472
473
474
475
476
477
478
	if ( input_file != NULL ) {

		full = read_reflections(input_file);
		if ( full == NULL ) {
			ERROR("Failed to read reflections from '%s'\n",
			      input_file);
			return 1;
		}
		free(input_file);
		if ( check_list_symmetry(full, sym) ) {
			ERROR("The input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
479
			      " have symmetry %s\n", symmetry_name(sym));
480
481
482
483
484
			return 1;
		}

	} else {
		random_intensities = 1;
485
486
	}

487
488
489
490
491
	if ( n < 1 ) {
		ERROR("Number of patterns must be at least 1.\n");
		return 1;
	}

492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
	if ( output_file == NULL ) {
		ERROR("You must pgive a filename for the output.\n");
		return 1;
	}
	ofh = fopen(output_file, "w");
	if ( ofh == NULL ) {
		ERROR("Couldn't open output file '%s'\n", output_file);
		return 1;
	}
	free(output_file);
	write_stream_header(ofh, argc, argv);

	image.det = det;
	image.width = det->max_fs;
	image.height = det->max_ss;

	image.lambda = ph_en_to_lambda(eV_to_J(beam->photon_energy));
	image.div = beam->divergence;
	image.bw = beam->bandwidth;
511
	image.profile_radius = 0.003e9;
512
	image.filename = malloc(256);
513
	image.copyme = NULL;
514

515
516
517
518
	if ( random_intensities ) {
		full = reflist_new();
	}

519
520
521
522
	qargs.full = full;
	pthread_mutex_init(&qargs.full_lock, NULL);
	qargs.n_to_do = n;
	qargs.n_done = 0;
523
	qargs.n_started = 0;
524
525
526
527
528
	qargs.sym = sym;
	qargs.random_intensities = random_intensities;
	qargs.cell = cell;
	qargs.template_image = &image;
	qargs.stream = ofh;
529
	qargs.cnoise = cnoise;
530
531
532
533
534
	qargs.max_q = largest_q(&image);

	for ( i=0; i<NBINS; i++ ) {
		qargs.n_ref[i] = 0;
		qargs.p_hist[i] = 0.0;
535
		qargs.p_max[i] = 0.0;
536
	}
537
538

	run_threads(n_threads, run_job, create_job, finalise_job,
539
	            &qargs, n, 0, 0, 0);
540

541
	if ( random_intensities ) {
542
		STATUS("Writing full intensities to %s\n", save_file);
543
		write_reflist(save_file, full);
544
545
	}

546
	if ( phist_file != NULL ) {
Thomas White's avatar
Thomas White committed
547

548
549
		fh = fopen(phist_file, "w");

Thomas White's avatar
Thomas White committed
550
		if ( fh != NULL ) {
551

Thomas White's avatar
Thomas White committed
552
553
			fprintf(fh, "1/d_nm #refl  pmean  pmax\n");

Thomas White's avatar
Thomas White committed
554
			for ( i=0; i<NBINS; i++ ) {
555

Thomas White's avatar
Thomas White committed
556
				double rcen;
557

Thomas White's avatar
Thomas White committed
558
559
560
561
562
563
564
565
				rcen = i/(double)NBINS*qargs.max_q
					  + qargs.max_q/(2.0*NBINS);
				fprintf(fh, "%.2f %7li %.3f %.3f\n", rcen/1.0e9,
					qargs.n_ref[i],
					qargs.p_hist[i]/qargs.n_ref[i],
					qargs.p_max[i]);

			}
566

Thomas White's avatar
Thomas White committed
567
568
569
570
571
572
			fclose(fh);

		} else {
			ERROR("Failed to open file '%s' for writing.\n",
			      phist_file);
		}
573
574
575

	}

576
577
578
579
	fclose(ofh);
	cell_free(cell);
	free_detector_geometry(det);
	free(beam);
Thomas White's avatar
Thomas White committed
580
	free_symoplist(sym);
581
	reflist_free(full);
582
	free(image.filename);
583
584
585

	return 0;
}