partial_sim.c 10.4 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
/*
 * partial_sim.c
 *
 * Generate partials for testing scaling
 *
 * (c) 2006-2011 Thomas White <taw@physics.org>
 *
 * 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 <assert.h>
24
#include <pthread.h>
25
26
27
28
29
30
31
32

#include "utils.h"
#include "reflist-utils.h"
#include "symmetry.h"
#include "beam-parameters.h"
#include "detector.h"
#include "geometry.h"
#include "stream.h"
33
#include "thread-pool.h"
34
35


36
static void mess_up_cell(UnitCell *cell, double cnoise)
37
38
39
40
41
{
	double ax, ay, az;
	double bx, by, bz;
	double cx, cy, cz;

Thomas White's avatar
Thomas White committed
42
43
	//STATUS("Real:\n");
	//cell_print(cell);
Thomas White's avatar
Thomas White committed
44

45
	cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
46
47
48
49
50
51
52
53
54
	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);
55
	cell_set_reciprocal(cell, ax, ay, az, bx, by, bz, cx, cy, cz);
Thomas White's avatar
Thomas White committed
56

Thomas White's avatar
Thomas White committed
57
58
	//STATUS("Changed:\n");
	//cell_print(cell);
59
60
}

61

62
63
/* For each reflection in "partial", fill in what the intensity would be
 * according to "full" */
64
static void calculate_partials(RefList *partial, double osf,
Thomas White's avatar
Thomas White committed
65
                               RefList *full, const SymOpList *sym,
66
67
                               int random_intensities,
                               pthread_mutex_t *full_lock)
68
69
70
71
72
73
{
	Reflection *refl;
	RefListIterator *iter;

	for ( refl = first_refl(partial, &iter);
	      refl != NULL;
Thomas White's avatar
Thomas White committed
74
75
	      refl = next_refl(refl, iter) )
	{
76
77
		signed int h, k, l;
		Reflection *rfull;
Thomas White's avatar
Thomas White committed
78
		double p, Ip, If;
79
80

		get_indices(refl, &h, &k, &l);
Thomas White's avatar
Thomas White committed
81
		get_asymm(sym, h, k, l, &h, &k, &l);
82
83
		p = get_partiality(refl);

84
		pthread_mutex_lock(full_lock);
85
		rfull = find_refl(full, h, k, l);
86
87
		pthread_mutex_unlock(full_lock);

88
		if ( rfull == NULL ) {
89
			if ( random_intensities ) {
90
91
92
93
94
95

				/* 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);
96
97
98
99
				rfull = add_refl(full, h, k, l);
				If = fabs(gaussian_noise(0.0, 1000.0));
				set_int(rfull, If);
				set_redundancy(rfull, 1);
100
101
				pthread_mutex_unlock(full_lock);

102
103
104
105
			} else {
				set_redundancy(refl, 0);
				If = 0.0;
			}
106
107
		} else {
			If = get_intensity(rfull);
108
109
110
111
			if ( random_intensities ) {
				int red = get_redundancy(rfull);
				set_redundancy(rfull, red+1);
			}
112
		}
113
114

		Ip = osf * p * If;
115

116
		Ip = gaussian_noise(Ip, 100.0);
117

118
		set_int(refl, Ip);
119
		set_esd_intensity(refl, 100.0);
120
121
122
123
124
125
126
127
128
129
130
131
132
	}
}


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"
133
134
135
136
137
138
" -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"
139
"\n"
140
141
142
" -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"
143
144
145
146
" -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"
"\n"
147
148
149
150
);
}


151
152
153
154
155
156
157
158
159
160
161
struct queue_args
{
	RefList *full;
	pthread_mutex_t full_lock;

	int n_done;
	int n_to_do;

	SymOpList *sym;
	int random_intensities;
	UnitCell *cell;
162
	double cnoise;
163
164
165
166
167
168
169
170
171
172
173
174
175
176
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

	struct image *template_image;

	FILE *stream;
};


struct worker_args
{
	struct queue_args *qargs;
	struct image image;
};


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

	wargs = malloc(sizeof(struct worker_args));

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

	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;

	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);
	calculate_partials(wargs->image.reflections, osf, qargs->full,
	                   qargs->sym, qargs->random_intensities,
	                   &qargs->full_lock);

	/* Give a slightly incorrect cell in the stream */
212
	mess_up_cell(wargs->image.indexed_cell, qargs->cnoise);
213
214
215
216
217
218
219
220
}


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

221
	write_chunk(qargs->stream, &wargs->image, NULL, STREAM_INTEGRATED);
222
223
224
225
226
227
228
229
230
231

	reflist_free(wargs->image.reflections);
	cell_free(wargs->image.indexed_cell);
	free(wargs);

	qargs->n_done++;
	progress_bar(qargs->n_done, qargs->n_to_do, "Simulating");
}


232
233
234
235
236
237
238
239
240
241
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;
242
	RefList *full = NULL;
Thomas White's avatar
Thomas White committed
243
244
	char *sym_str = NULL;
	SymOpList *sym;
245
246
	UnitCell *cell = NULL;
	FILE *ofh;
247
	int n = 2;
248
	int random_intensities = 0;
249
	char *save_file = NULL;
250
251
252
	struct queue_args qargs;
	struct image image;
	int n_threads = 1;
253
254
	double cnoise = 0.0;
	char *rval;
255
256
257
258
259
260
261
262
263
264

	/* 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'},
265
		{"save-random",        1, NULL,               'r'},
266
		{"cnoise",             1, NULL,               'c'},
267
268
269
270
		{0, 0, NULL, 0}
	};

	/* Short options */
271
	while ((c = getopt_long(argc, argv, "hi:o:b:p:g:y:n:r:j:c:",
Thomas White's avatar
Thomas White committed
272
273
	                        longopts, NULL)) != -1)
	{
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
		switch (c) {
		case 'h' :
			show_help(argv[0]);
			return 0;

		case 'o' :
			output_file = strdup(optarg);
			break;

		case 'i' :
			input_file = strdup(optarg);
			break;

		case 'b' :
			beamfile = strdup(optarg);
			break;

		case 'p' :
			cellfile = strdup(optarg);
			break;

		case 'g' :
			geomfile = strdup(optarg);
			break;

		case 'y' :
Thomas White's avatar
Thomas White committed
300
			sym_str = strdup(optarg);
301
302
			break;

303
304
305
306
		case 'n' :
			n = atoi(optarg);
			break;

307
308
309
310
		case 'r' :
			save_file = strdup(optarg);
			break;

311
312
313
314
		case 'j' :
			n_threads = atoi(optarg);
			break;

315
316
317
318
319
320
321
322
		case 'c' :
			cnoise = strtod(optarg, &rval);
			if ( *rval != '\0' ) {
				ERROR("Invalid cell noise value.\n");
				return 1;
			}
			break;

323
324
325
326
327
328
329
330
		case 0 :
			break;

		default :
			return 1;
		}
	}

331
332
333
334
335
	if ( n_threads < 1 ) {
		ERROR("Invalid number of threads.\n");
		return 1;
	}

336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
	/* 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
360
361
362
363
364
365
	if ( !cell_is_sensible(cell) ) {
		ERROR("Invalid unit cell parameters:\n");
		cell_print(cell);
		return 1;
	}

366
367
368
369
370
371
	/* Load geometry */
	if ( geomfile == NULL ) {
		ERROR("You need to give a geometry file.\n");
		return 1;
	}
	det = get_detector_geometry(geomfile);
372
	if ( det == NULL ) {
373
374
375
376
377
		ERROR("Failed to read geometry from '%s'\n", geomfile);
		return 1;
	}
	free(geomfile);

Thomas White's avatar
Thomas White committed
378
379
380
381
	if ( sym_str == NULL ) sym_str = strdup("1");
	sym = get_pointgroup(sym_str);
	free(sym_str);

382
383
	if ( save_file == NULL ) save_file = strdup("partial_sim.hkl");

384
	/* Load (full) reflections */
385
386
387
388
389
390
391
392
393
394
395
	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
396
			      " have symmetry %s\n", symmetry_name(sym));
397
398
399
400
401
			return 1;
		}

	} else {
		random_intensities = 1;
402
403
	}

404
405
406
407
408
	if ( n < 1 ) {
		ERROR("Number of patterns must be at least 1.\n");
		return 1;
	}

409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
	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;
Thomas White's avatar
Thomas White committed
428
	image.profile_radius = 0.003e9;
429
	image.i0_available = 0;
430
431
	image.filename = malloc(256);

432
433
434
435
	if ( random_intensities ) {
		full = reflist_new();
	}

436
437
438
439
440
441
442
443
444
	qargs.full = full;
	pthread_mutex_init(&qargs.full_lock, NULL);
	qargs.n_to_do = n;
	qargs.n_done = 0;
	qargs.sym = sym;
	qargs.random_intensities = random_intensities;
	qargs.cell = cell;
	qargs.template_image = &image;
	qargs.stream = ofh;
445
	qargs.cnoise = cnoise;
446
447

	run_threads(n_threads, run_job, create_job, finalise_job,
448
	            &qargs, n, 0, 0, 0);
449

450
	if ( random_intensities ) {
451
452
		STATUS("Writing full intensities to %s\n", save_file);
		write_reflist(save_file, full, cell);
453
454
	}

455
456
457
458
	fclose(ofh);
	cell_free(cell);
	free_detector_geometry(det);
	free(beam);
Thomas White's avatar
Thomas White committed
459
	free_symoplist(sym);
460
	reflist_free(full);
461
	free(image.filename);
462
463
464

	return 0;
}