compare_hkl.c 17.8 KB
Newer Older
1
2
3
4
5
/*
 * compare_hkl.c
 *
 * Compare reflection lists
 *
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>
41
#include <gsl/gsl_errno.h>
Thomas White's avatar
Thomas White committed
42
#include <gsl/gsl_statistics.h>
43
44

#include "utils.h"
Thomas White's avatar
Thomas White committed
45
#include "statistics.h"
46
#include "symmetry.h"
Thomas White's avatar
Thomas White committed
47
#include "reflist-utils.h"
48
#include "cell-utils.h"
49
50


Thomas White's avatar
Thomas White committed
51
enum fom
52
{
Thomas White's avatar
Thomas White committed
53
54
55
56
57
58
	FOM_R1I,
	FOM_R1F,
	FOM_R2,
	FOM_RSPLIT,
	FOM_CC,
	FOM_CCSTAR
59
60
61
};


Thomas White's avatar
Thomas White committed
62
static enum fom get_fom(const char *s)
63
{
Thomas White's avatar
Thomas White committed
64
65
66
67
68
69
70
71
	if ( strcasecmp(s, "r1i") == 0 ) return FOM_R1I;
	if ( strcasecmp(s, "r1f") == 0 ) return FOM_R1F;
	if ( strcasecmp(s, "r2") == 0 ) return FOM_R2;
	if ( strcasecmp(s, "rsplit") == 0 ) return FOM_RSPLIT;
	if ( strcasecmp(s, "cc") == 0 ) return FOM_CC;
	if ( strcasecmp(s, "ccstar") == 0 ) return FOM_CCSTAR;

	ERROR("Unknown figure of merit '%s'.\n", s);
72
73
74
75
	exit(1);
}


76
77
static void show_help(const char *s)
{
78
	printf("Syntax: %s [options] <file1.hkl> <file2.hkl>\n\n", s);
79
80
81
	printf(
"Compare intensity lists.\n"
"\n"
82
"  -y, --symmetry=<sym>       The symmetry of both the input files.\n"
Thomas White's avatar
Thomas White committed
83
"  -p, --pdb=<filename>       PDB file to use.\n"
Thomas White's avatar
Thomas White committed
84
85
86
87
88
"      --fom=<FoM>            Calculate this figure of merit  Choose from:.\n"
"                              R1I, R1F, R2, Rsplit,\n"
"                              CCF, CCI, CCFstar, CCIstar.\n"
"      --nshells=<n>          Use <n> resolution shells.\n"
"  -u                         Force scale factor to 1.\n"
89
"      --shell-file=<file>    Write resolution shells to <file>.\n"
Thomas White's avatar
Thomas White committed
90
91
92
93
94
"\n"
"You can control which reflections are included in the calculation:\n"
"\n"
"      --ignore-negs          Ignore reflections with negative intensities.\n"
"      --zero-negs            Set negative intensities to zero.\n"
95
"      --sigma-cutoff=<n>     Discard reflections with I/sigma(I) < n.\n"
Thomas White's avatar
Thomas White committed
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
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
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
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
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
"      --rmin=<res>           Set a lower resolution limit (m^-1).\n"
"      --rmax=<res>           Set an upper resolution limit (m^-1).\n"
"\n"
"  -h, --help                 Display this help message.\n"
);
}


struct fom_context
{
	enum fom fom;
	int nshells;

	/* For R-factors */
	double *num;
	double *den;

	/* For CCs */
	double **vec1;
	double **vec2;
	int *n;
	int nmax;
};


static struct fom_context *init_fom(enum fom fom, int nmax, int nshells)
{
	struct fom_context *fctx;
	int i;

	fctx = malloc(sizeof(struct fom_context));
	if ( fctx == NULL ) return NULL;

	fctx->fom = fom;
	fctx->nshells = nshells;

	switch ( fctx->fom ) {

		case FOM_R1I :
		case FOM_R1F :
		case FOM_R2 :
		case FOM_RSPLIT :
		fctx->num = malloc(nshells*sizeof(double));
		fctx->den = malloc(nshells*sizeof(double));
		if ( (fctx->num == NULL) || (fctx->den == NULL) ) return NULL;
		for ( i=0; i<nshells; i++ ) {
			fctx->num[i] = 0.0;
			fctx->den[i] = 0.0;
		}
		break;

		case FOM_CC :
		case FOM_CCSTAR :
		fctx->vec1 = malloc(nshells*sizeof(double *));
		fctx->vec2 = malloc(nshells*sizeof(double *));
		if ( (fctx->vec1 == NULL) || (fctx->vec2 == NULL) ) return NULL;
		for ( i=0; i<nshells; i++ ) {
			fctx->vec1[i] = malloc(nmax*sizeof(double));
			if ( fctx->vec1[i] == NULL ) return NULL;
			fctx->vec2[i] = malloc(nmax*sizeof(double));
			if ( fctx->vec2[i] == NULL ) return NULL;
			fctx->n = malloc(nshells*sizeof(int));
			if ( fctx->n == NULL ) return NULL;
		}
		for ( i=0; i<nshells; i++ ) {
			fctx->n[i] = 0;
		}

		fctx->nmax = nmax;
		break;

	}

	return fctx;
}


static void add_to_fom(struct fom_context *fctx, double i1, double i2,
                       int bin)
{
	double f1, f2;

	/* Negative intensities have already been weeded out. */
	f1 = sqrt(i1);
	f2 = sqrt(i2);

	switch ( fctx->fom ) {

		case FOM_R1I :
		fctx->num[bin] += fabs(i1 - i2);
		fctx->den[bin] += i1;
		break;

		case FOM_R1F :
		fctx->num[bin] += fabs(f1 - f2);
		fctx->den[bin] += f1;
		break;

		case FOM_R2 :
		fctx->num[bin] += pow(i1 - i2, 2.0);
		fctx->den[bin] += pow(i1, 2.0);
		break;

		case FOM_RSPLIT :
		fctx->num[bin] += fabs(i1 - i2);
		fctx->den[bin] += i1 + i2;
		break;

		case FOM_CC :
		case FOM_CCSTAR :
		assert(fctx->n[bin] < fctx->nmax);
		fctx->vec1[bin][fctx->n[bin]] = i1;
		fctx->vec2[bin][fctx->n[bin]] = i2;
		fctx->n[bin]++;
		break;

	}
}


static double fom_overall(struct fom_context *fctx)
{
	double overall_num = INFINITY;
	double overall_den = 0.0;
	int i;
	double *overall_vec1;
	double *overall_vec2;
	int overall_n;
	double cc = INFINITY;

	switch ( fctx->fom ) {

		case FOM_R1I :
		case FOM_R1F :
		case FOM_R2 :
		case FOM_RSPLIT :
		overall_num = 0.0;
		overall_den = 0.0;
		for ( i=0; i<fctx->nshells; i++ ) {
			overall_num += fctx->num[i];
			overall_den += fctx->den[i];
		}
		break;

		case FOM_CC :
		case FOM_CCSTAR :
		overall_vec1 = malloc(fctx->nmax*sizeof(double));
		overall_vec2 = malloc(fctx->nmax*sizeof(double));
		overall_n = 0;
		for ( i=0; i<fctx->nshells; i++ ) {
			int j;
			for ( j=0; j<fctx->n[i]; j++ ) {
				overall_vec1[overall_n] = fctx->vec1[i][j];
				overall_vec2[overall_n] = fctx->vec2[i][j];
				overall_n++;
			}
		}
		cc = gsl_stats_correlation(overall_vec1, 1, overall_vec2, 1,
		                           overall_n);
		free(overall_vec1);
		free(overall_vec2);
		break;

	}

	switch ( fctx->fom ) {

		case FOM_R1I :
		case FOM_R1F :
		return overall_num/overall_den;

		case FOM_R2 :
		return sqrt(overall_num/overall_den);

		case FOM_RSPLIT :
		return 2.0*(overall_num/overall_den) / sqrt(2.0);

		case FOM_CC :
		return cc;

		case FOM_CCSTAR :
		return sqrt((2.0*cc)/(1.0+cc));

	}

	ERROR("This point is never reached.\n");
	abort();
}


static double fom_shell(struct fom_context *fctx, int i)
{
	double cc;

	switch ( fctx->fom ) {

		case FOM_R1I :
		case FOM_R1F :
		return fctx->num[i]/fctx->den[i];

		case FOM_R2 :
		return sqrt(fctx->num[i]/fctx->den[i]);

		case FOM_RSPLIT :
		return 2.0*(fctx->num[i]/fctx->den[i]) / sqrt(2.0);

		case FOM_CC :
		return gsl_stats_correlation(fctx->vec1[i], 1, fctx->vec2[i], 1,
		                             fctx->n[i]);

		case FOM_CCSTAR :
		cc = gsl_stats_correlation(fctx->vec1[i], 1, fctx->vec2[i], 1,
		                           fctx->n[i]);
		return sqrt((2.0*cc)/(1.0+cc));

	}

	ERROR("This point is never reached.\n");
	abort();
315
316
317
}


Thomas White's avatar
Thomas White committed
318
319
static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
                   double rmin, double rmax, enum fom fom,
320
                   int config_unity, int nshells, const char *filename)
Thomas White's avatar
Thomas White committed
321
{
Thomas White's avatar
Thomas White committed
322
323
324
	int *cts;
	double *rmins;
	double *rmaxs;
325
	double total_vol, vol_per_shell;
Thomas White's avatar
Thomas White committed
326
	int i;
Thomas White's avatar
Thomas White committed
327
328
	Reflection *refl1;
	RefListIterator *iter;
329
	FILE *fh;
Thomas White's avatar
Thomas White committed
330
331
332
333
334
335
336
	double scale;
	struct fom_context *fctx;

	cts = malloc(nshells*sizeof(int));
	rmins = malloc(nshells*sizeof(double));
	rmaxs = malloc(nshells*sizeof(double));
	fctx = init_fom(fom, num_reflections(list1), nshells);
Thomas White's avatar
Thomas White committed
337

Thomas White's avatar
Thomas White committed
338
339
340
	if ( (fctx==NULL) || (cts==NULL) || (rmins==NULL) || (rmaxs==NULL) )
	{
		ERROR("Couldn't allocate memory for resolution shells.\n");
Thomas White's avatar
Thomas White committed
341
342
343
		return;
	}

Thomas White's avatar
Thomas White committed
344
	for ( i=0; i<nshells; i++ ) {
Thomas White's avatar
Thomas White committed
345
		cts[i] = 0;
Thomas White's avatar
Thomas White committed
346
347
	}

Thomas White's avatar
Thomas White committed
348
349
350
351
352
	if ( config_unity ) {
		scale = 1.0;
	} else {
		scale = stat_scale_intensity(list1, list2);
	}
353

Thomas White's avatar
Thomas White committed
354
	/* Calculate the resolution bins */
355
	total_vol = pow(rmax, 3.0) - pow(rmin, 3.0);
Thomas White's avatar
Thomas White committed
356
	vol_per_shell = total_vol / nshells;
357
	rmins[0] = rmin;
Thomas White's avatar
Thomas White committed
358
	for ( i=1; i<nshells; i++ ) {
359
360
361
362
363
364

		double r;

		r = vol_per_shell + pow(rmins[i-1], 3.0);
		r = pow(r, 1.0/3.0);

365
		/* Shells of constant volume */
366
367
368
		rmaxs[i-1] = r;
		rmins[i] = r;

369
370
371
372
		/* Shells of constant thickness */
		//rmins[i] = rmins[i-1] + (rmax-rmin)/NBINS;
		//rmaxs[i-1] = rmins[i-1] + (rmax-rmin)/NBINS;

373
	}
Thomas White's avatar
Thomas White committed
374
	rmaxs[nshells-1] = rmax;
375

Thomas White's avatar
Thomas White committed
376
377
	for ( refl1 = first_refl(list1, &iter);
	      refl1 != NULL;
378
379
	      refl1 = next_refl(refl1, iter) )
	{
Thomas White's avatar
Thomas White committed
380
		signed int h, k, l;
381
		double d;
Thomas White's avatar
Thomas White committed
382
		int bin;
Thomas White's avatar
Thomas White committed
383
		double i1, i2;
384
		int j;
385
		Reflection *refl2;
Thomas White's avatar
Thomas White committed
386

Thomas White's avatar
Thomas White committed
387
		get_indices(refl1, &h, &k, &l);
388
		d = 2.0 * resolution(cell, h, k, l);
Thomas White's avatar
Thomas White committed
389

390
		bin = -1;
Thomas White's avatar
Thomas White committed
391
		for ( j=0; j<nshells; j++ ) {
392
393
394
395
396
			if ( (d>rmins[j]) && (d<=rmaxs[j]) ) {
				bin = j;
				break;
			}
		}
397

Thomas White's avatar
Thomas White committed
398
399
400
		/* Allow for slight rounding errors */
		if ( (bin == -1) && (d <= rmins[0]) ) bin = 0;
		assert(bin != -1);
Thomas White's avatar
Thomas White committed
401

402
403
404
		refl2 = find_refl(list2, h, k, l);
		if ( refl2 == NULL ) continue;

Thomas White's avatar
Thomas White committed
405
		i1 = get_intensity(refl1);
Thomas White's avatar
Thomas White committed
406
		i2 = scale * get_intensity(refl2);
407

Thomas White's avatar
Thomas White committed
408
409
410
		add_to_fom(fctx, i1, i2, bin);
		cts[bin]++;
	}
411

Thomas White's avatar
Thomas White committed
412
	switch ( fom ) {
413

Thomas White's avatar
Thomas White committed
414
415
416
		case FOM_R1I :
		STATUS("Overall R1(I) = %.2f %%\n", 100.0*fom_overall(fctx));
		break;
417

Thomas White's avatar
Thomas White committed
418
419
420
		case FOM_R1F :
		STATUS("Overall R1(F) = %.2f %%\n", 100.0*fom_overall(fctx));
		break;
421

Thomas White's avatar
Thomas White committed
422
423
424
		case FOM_R2 :
		STATUS("Overall R(2) = %.2f %%\n", 100.0*fom_overall(fctx));
		break;
425

Thomas White's avatar
Thomas White committed
426
427
428
		case FOM_RSPLIT :
		STATUS("Overall Rsplit = %.2f %%\n", 100.0*fom_overall(fctx));
		break;
Thomas White's avatar
Thomas White committed
429

Thomas White's avatar
Thomas White committed
430
431
432
433
434
435
436
		case FOM_CC :
		STATUS("Overall CC = %.7f\n", fom_overall(fctx));
		break;

		case FOM_CCSTAR :
		STATUS("Overall CC* = %.7f\n", fom_overall(fctx));
		break;
Thomas White's avatar
Thomas White committed
437

438
439
	}

440
	fh = fopen(filename, "w");
441
	if ( fh == NULL ) {
442
		ERROR("Couldn't open '%s'\n", filename);
443
444
445
		return;
	}

Thomas White's avatar
Thomas White committed
446
447
448
	switch ( fom ) {

		case FOM_R1I :
Thomas White's avatar
Thomas White committed
449
		fprintf(fh, "1/d centre    R1(I)/%%       nref      d / A\n");
Thomas White's avatar
Thomas White committed
450
		break;
451

Thomas White's avatar
Thomas White committed
452
		case FOM_R1F :
Thomas White's avatar
Thomas White committed
453
		fprintf(fh, "1/d centre    R1(F)/%%       nref      d / A\n");
454
455
		break;

Thomas White's avatar
Thomas White committed
456
		case FOM_R2 :
Thomas White's avatar
Thomas White committed
457
		fprintf(fh, "1/d centre       R2/%%       nref      d / A\n");
458
459
		break;

Thomas White's avatar
Thomas White committed
460
		case FOM_RSPLIT :
Thomas White's avatar
Thomas White committed
461
		fprintf(fh, "1/d centre   Rsplit/%%       nref      d / A\n");
462
463
		break;

Thomas White's avatar
Thomas White committed
464
		case FOM_CC :
Thomas White's avatar
Thomas White committed
465
		fprintf(fh, "1/d centre         CC       nref      d / A\n");
Thomas White's avatar
Thomas White committed
466
467
468
		break;

		case FOM_CCSTAR :
Thomas White's avatar
Thomas White committed
469
		fprintf(fh, "1/d centre        CC*       nref      d / A\n");
470
471
472
		break;

	}
473

Thomas White's avatar
Thomas White committed
474
	for ( i=0; i<nshells; i++ ) {
Thomas White's avatar
Thomas White committed
475

476
		double r, cen;
477
		cen = rmins[i] + (rmaxs[i] - rmins[i])/2.0;
478

Thomas White's avatar
Thomas White committed
479
		r = fom_shell(fctx, i);
480

Thomas White's avatar
Thomas White committed
481
		switch ( fom ) {
482

Thomas White's avatar
Thomas White committed
483
484
485
486
487
488
		case FOM_R1I :
		case FOM_R1F :
		case FOM_R2 :
		case FOM_RSPLIT :
		fprintf(fh, "%10.3f %10.2f %10i %10.2f\n",
		        cen*1.0e-9, r*100.0, cts[i], (1.0/cen)*1e10);
489

Thomas White's avatar
Thomas White committed
490
		break;
491

Thomas White's avatar
Thomas White committed
492
493
494
495
		case FOM_CC :
		case FOM_CCSTAR :
		fprintf(fh, "%10.3f %10.7f %10i %10.2f\n",
		        cen*1.0e-9, r, cts[i], (1.0/cen)*1e10);
496

Thomas White's avatar
Thomas White committed
497
		break;
Thomas White's avatar
Thomas White committed
498
499
500

	}

Thomas White's avatar
Thomas White committed
501
	}
Thomas White's avatar
Thomas White committed
502

Thomas White's avatar
Thomas White committed
503
	fclose(fh);
Thomas White's avatar
Thomas White committed
504
505
506
}


507
508
509
int main(int argc, char *argv[])
{
	int c;
510
	UnitCell *cell;
511
512
	char *afile = NULL;
	char *bfile = NULL;
Thomas White's avatar
Thomas White committed
513
514
	char *sym_str = NULL;
	SymOpList *sym;
Thomas White's avatar
Thomas White committed
515
516
517
	int ncom, nrej, nneg, nres;
	RefList *list1_acc;
	RefList *list2_acc;
Thomas White's avatar
Thomas White committed
518
519
	RefList *list1;
	RefList *list2;
520
521
	RefList *list1_raw;
	RefList *list2_raw;
Thomas White's avatar
Thomas White committed
522
	enum fom fom = FOM_R1I;
523
	char *pdb = NULL;
524
525
	float rmin_fix = -1.0;
	float rmax_fix = -1.0;
Thomas White's avatar
Thomas White committed
526
	double rmin, rmax;
Thomas White's avatar
Thomas White committed
527
528
	Reflection *refl1;
	RefListIterator *iter;
529
	float sigma_cutoff = -INFINITY;
Thomas White's avatar
Thomas White committed
530
531
532
533
	int config_ignorenegs = 0;
	int config_zeronegs = 0;
	int config_unity = 0;
	int nshells = 10;
534
	char *shell_file = NULL;
535
536
537
538

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
539
		{"symmetry",           1, NULL,               'y'},
540
		{"pdb",                1, NULL,               'p'},
541
542
		{"rmin",               1, NULL,                2},
		{"rmax",               1, NULL,                3},
Thomas White's avatar
Thomas White committed
543
		{"fom",                1, NULL,                4},
544
		{"sigma-cutoff",       1, NULL,                5},
Thomas White's avatar
Thomas White committed
545
		{"nshells",            1, NULL,                6},
546
		{"shell-file",         1, NULL,                7},
Thomas White's avatar
Thomas White committed
547
548
		{"ignore-negs",        0, &config_ignorenegs,  1},
		{"zero-negs",          0, &config_zeronegs,    1},
549
550
551
552
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
553
	while ((c = getopt_long(argc, argv, "hy:p:u",
554
555
	                        longopts, NULL)) != -1)
	{
556
557

		switch (c) {
Thomas White's avatar
Thomas White committed
558
559

			case 'h' :
560
561
562
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
563
			case 'y' :
Thomas White's avatar
Thomas White committed
564
			sym_str = strdup(optarg);
565
566
			break;

Thomas White's avatar
Thomas White committed
567
			case 'p' :
568
569
570
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
571
			case 'u' :
572
573
574
			config_unity = 1;
			break;

Thomas White's avatar
Thomas White committed
575
			case 0 :
576
577
			break;

Thomas White's avatar
Thomas White committed
578
			case 2 :
579
580
581
582
583
584
			if ( sscanf(optarg, "%e", &rmin_fix) != 1 ) {
				ERROR("Invalid value for --rmin\n");
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
585
			case 3 :
586
587
588
589
590
591
			if ( sscanf(optarg, "%e", &rmax_fix) != 1 ) {
				ERROR("Invalid value for --rmax\n");
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
592
			case 4 :
Thomas White's avatar
Thomas White committed
593
			fom = get_fom(optarg);
594
595
			break;

596
597
598
599
600
601
602
			case 5 :
			if ( sscanf(optarg, "%f", &sigma_cutoff) != 1 ) {
				ERROR("Invalid value for --sigma-cutoff\n");
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
603
604
605
606
607
608
609
			case 6 :
			if ( sscanf(optarg, "%i", &nshells) != 1 ) {
				ERROR("Invalid value for --nshells\n");
				return 1;
			}
			break;

610
611
612
613
			case 7 :
			shell_file = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
614
			default :
615
616
			ERROR("Unhandled option '%c'\n", c);
			break;
Thomas White's avatar
Thomas White committed
617

618
619
620
621
		}

	}

622
623
624
625
626
	if ( argc != (optind+2) ) {
		ERROR("Please provide exactly two HKL files to compare.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
	if ( !config_ignorenegs && !config_zeronegs ) {
		switch ( fom )
		{
			case FOM_R1F :
			ERROR("Your chosen figure of merit involves converting"
			      " intensities to structure factors, but you have"
			      " not specified how to handle negative"
			      " intensities.\n");
			ERROR("Please try again with --ignore-negs or"
			      " --zero-negs.\n");
			exit(1);

			case FOM_R2 :
			case FOM_R1I :
			case FOM_RSPLIT :
			case FOM_CC :
			case FOM_CCSTAR :
			break;
		}
	}

Thomas White's avatar
Thomas White committed
648
649
	if ( sym_str == NULL ) {
		sym_str = strdup("1");
650
	}
Thomas White's avatar
Thomas White committed
651
652
	sym = get_pointgroup(sym_str);
	free(sym_str);
653

654
655
656
	afile = strdup(argv[optind++]);
	bfile = strdup(argv[optind]);

657
	if ( pdb == NULL ) {
Thomas White's avatar
Thomas White committed
658
659
		ERROR("You must provide a PDB file with the unit cell.\n");
		exit(1);
660
661
	}

662
663
	if ( shell_file == NULL ) shell_file = strdup("shells.dat");

664
665
666
	cell = load_cell_from_pdb(pdb);
	free(pdb);

667
668
	list1_raw = read_reflections(afile);
	if ( list1_raw == NULL ) {
Thomas White's avatar
Thomas White committed
669
		ERROR("Couldn't read file '%s'\n", afile);
670
671
		return 1;
	}
Thomas White's avatar
Thomas White committed
672

673
674
	list2_raw = read_reflections(bfile);
	if ( list2_raw == NULL ) {
Thomas White's avatar
Thomas White committed
675
		ERROR("Couldn't read file '%s'\n", bfile);
676
677
678
		return 1;
	}

Thomas White's avatar
Thomas White committed
679
	/* Check that the intensities have the correct symmetry */
680
	if ( check_list_symmetry(list1_raw, sym) ) {
Thomas White's avatar
Thomas White committed
681
		ERROR("The first input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
682
		      " have symmetry %s\n", symmetry_name(sym));
Thomas White's avatar
Thomas White committed
683
684
		return 1;
	}
685
	if ( check_list_symmetry(list2_raw, sym) ) {
Thomas White's avatar
Thomas White committed
686
		ERROR("The second input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
687
		      " have symmetry %s\n", symmetry_name(sym));
Thomas White's avatar
Thomas White committed
688
689
690
		return 1;
	}

Thomas White's avatar
Thomas White committed
691
692
693
694
695
696
697
698
699
700
701
702
	resolution_limits(list1_raw, cell, &rmin, &rmax);
	STATUS("%s: %i reflections, resolution range %.2f to %.2f Angstroms"
	       " (%.5f to %.5f nm^-1).\n", afile,
	       num_reflections(list1_raw),
	       1e10/rmin, 1e10/rmax, rmin/1e9, rmax/1e9);

	resolution_limits(list2_raw, cell, &rmin, &rmax);
	STATUS("%s: %i reflections, resolution range %.2f to %.2f Angstroms"
	       " (%.5f to %.5f nm^-1).\n", bfile,
	       num_reflections(list2_raw),
	       1e10/rmin, 1e10/rmax, rmin/1e9, rmax/1e9);

703
704
705
	list1 = asymmetric_indices(list1_raw, sym);
	list2 = asymmetric_indices(list2_raw, sym);

Thomas White's avatar
Thomas White committed
706
	/* Select reflections to be used */
707
	ncom = 0;
708
	nrej = 0;
Thomas White's avatar
Thomas White committed
709
710
711
712
	nneg = 0;
	nres = 0;
	list1_acc = reflist_new();
	list2_acc = reflist_new();
Thomas White's avatar
Thomas White committed
713
714
	for ( refl1 = first_refl(list1, &iter);
	      refl1 != NULL;
Thomas White's avatar
Thomas White committed
715
716
	      refl1 = next_refl(refl1, iter) )
	{
Thomas White's avatar
Thomas White committed
717
		signed int h, k, l;
718
		double val1, val2;
719
		double esd1, esd2;
Thomas White's avatar
Thomas White committed
720
		Reflection *refl2;
Thomas White's avatar
Thomas White committed
721
722
		Reflection *refl1_acc;
		Reflection *refl2_acc;
723

Thomas White's avatar
Thomas White committed
724
		get_indices(refl1, &h, &k, &l);
725

726
727
		refl2 = find_refl(list2, h, k, l);
		if ( refl2 == NULL ) continue;
728

Thomas White's avatar
Thomas White committed
729
730
731
		val1 = get_intensity(refl1);
		val2 = get_intensity(refl2);

732
733
734
735
736
737
738
739
740
741
		esd1 = get_esd_intensity(refl1);
		esd2 = get_esd_intensity(refl2);

		if ( (val1 < sigma_cutoff * esd1)
		  || (val2 < sigma_cutoff * esd2) )
		{
			nrej++;
			continue;
		}

Thomas White's avatar
Thomas White committed
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
		if ( config_ignorenegs && ((val1 < 0.0) || (val2 < 0.0)) ) {
			nneg++;
			continue;
		}

		if ( config_zeronegs ) {
			int d = 0;
			if ( val1 < 0.0 ) {
				val1 = 0.0;
				d = 1;
			}
			if ( val2 < 0.0 ) {
				val2 = 0.0;
				d = 1;
			}
			if ( d ) nneg++;
		}

		if ( rmin_fix > 0.0 ) {
			double res = 2.0*resolution(cell, h, k, l);
			if ( res < rmin_fix ) {
				nres++;
				continue;
			}
		}

		if ( rmax_fix > 0.0 ) {
			double res = 2.0*resolution(cell, h, k, l);
Thomas White's avatar
Thomas White committed
770
			if ( res > rmax_fix ) {
Thomas White's avatar
Thomas White committed
771
772
773
774
775
776
777
778
779
780
781
782
783
				nres++;
				continue;
			}
		}

		refl1_acc = add_refl(list1_acc, h, k, l);
		copy_data(refl1_acc, refl1);
		set_intensity(refl1_acc, val1);

		refl2_acc = add_refl(list2_acc, h, k, l);
		copy_data(refl2_acc, refl2);
		set_intensity(refl2_acc, val2);

784
785
		ncom++;

786
	}
Thomas White's avatar
Thomas White committed
787

Thomas White's avatar
Thomas White committed
788
789
790
791
792
	gsl_set_error_handler_off();

	if ( nrej > 0 ) {
		STATUS("Discarded %i reflection pairs because either or both"
		       " versions had I/sigma(I) < %f.\n", nrej, sigma_cutoff);
793
794
	}

Thomas White's avatar
Thomas White committed
795
796
	if ( config_ignorenegs && (nneg > 0) ) {
		STATUS("Discarded %i reflection pairs because either or both"
Thomas White's avatar
Thomas White committed
797
		       " versions had negative intensities.\n", nneg);
Thomas White's avatar
Thomas White committed
798
	}
799

Thomas White's avatar
Thomas White committed
800
801
	if ( config_zeronegs && (nneg > 0) ) {
		STATUS("For %i reflection pairs, either or both versions had"
Thomas White's avatar
Thomas White committed
802
		       " negative intensities which were set to zero.\n", nneg);
803
804
	}

Thomas White's avatar
Thomas White committed
805
806
807
	if ( nres > 0 ) {
		STATUS("%i reflection pairs rejected because either or both "
		       " versions were outside the resolution range.\n", nres);
Thomas White's avatar
Thomas White committed
808
809
	}

Thomas White's avatar
Thomas White committed
810
811
812
813
814
815
816
817
818
819
820
821
822
823
	STATUS("%i reflection pairs accepted.\n", ncom);

	resolution_limits(list1_acc, cell, &rmin, &rmax);
	resolution_limits(list2_acc, cell, &rmin, &rmax);
	STATUS("Accepted resolution range: %f to %f nm^-1"
	       " (%.2f to %.2f Angstroms).\n",
	       rmin/1e9, rmax/1e9, 1e10/rmin, 1e10/rmax);

	reflist_free(list1_raw);
	reflist_free(list2_raw);
	reflist_free(list1);
	reflist_free(list2);

	do_fom(list1_acc, list2_acc, cell, rmin, rmax, fom, config_unity,
824
	       nshells, shell_file);
Thomas White's avatar
Thomas White committed
825

826
	free(shell_file);
Thomas White's avatar
Thomas White committed
827
828
829
	reflist_free(list1_acc);
	reflist_free(list2_acc);

830
831
	return 0;
}