compare_hkl.c 17.5 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
89
90
91
92
93
"      --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"
"\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"
94
"      --sigma-cutoff=<n>     Discard reflections with I/sigma(I) < n.\n"
Thomas White's avatar
Thomas White committed
95
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
"      --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();
314
315
316
}


Thomas White's avatar
Thomas White committed
317
318
319
static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
                   double rmin, double rmax, enum fom fom,
                   int config_unity, int nshells)
Thomas White's avatar
Thomas White committed
320
{
Thomas White's avatar
Thomas White committed
321
322
323
	int *cts;
	double *rmins;
	double *rmaxs;
324
	double total_vol, vol_per_shell;
Thomas White's avatar
Thomas White committed
325
	int i;
Thomas White's avatar
Thomas White committed
326
327
	Reflection *refl1;
	RefListIterator *iter;
328
	FILE *fh;
Thomas White's avatar
Thomas White committed
329
330
331
332
333
334
335
	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
336

Thomas White's avatar
Thomas White committed
337
338
339
	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
340
341
342
		return;
	}

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

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

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

		double r;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Thomas White's avatar
Thomas White committed
429
430
431
432
433
434
435
		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
436

437
438
	}

439
440
441
442
443
444
	fh = fopen("shells.dat", "w");
	if ( fh == NULL ) {
		ERROR("Couldn't open 'shells.dat'\n");
		return;
	}

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

		case FOM_R1I :
		fprintf(fh, "1/d centre  R1(I) / %%       nref      d (A)\n");
		break;
450

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

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

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

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

		case FOM_CCSTAR :
		fprintf(fh, "1/d centre        CC*       nref      d (A)\n");
469
470
471
		break;

	}
472

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

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

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

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

Thomas White's avatar
Thomas White committed
482
483
484
485
486
487
		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);
488

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

Thomas White's avatar
Thomas White committed
491
492
493
494
		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);
495

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

	}

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

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


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

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

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

		switch (c) {
Thomas White's avatar
Thomas White committed
555
556

			case 'h' :
557
558
559
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
560
			case 'y' :
Thomas White's avatar
Thomas White committed
561
			sym_str = strdup(optarg);
562
563
			break;

Thomas White's avatar
Thomas White committed
564
			case 'p' :
565
566
567
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
568
			case 'u' :
569
570
571
			config_unity = 1;
			break;

Thomas White's avatar
Thomas White committed
572
			case 0 :
573
574
			break;

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

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

Thomas White's avatar
Thomas White committed
589
			case 4 :
Thomas White's avatar
Thomas White committed
590
			fom = get_fom(optarg);
591
592
			break;

593
594
595
596
597
598
599
			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
600
601
602
603
604
605
606
			case 6 :
			if ( sscanf(optarg, "%i", &nshells) != 1 ) {
				ERROR("Invalid value for --nshells\n");
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
607
			default :
608
609
			ERROR("Unhandled option '%c'\n", c);
			break;
Thomas White's avatar
Thomas White committed
610

611
612
613
614
		}

	}

615
616
617
618
619
	if ( argc != (optind+2) ) {
		ERROR("Please provide exactly two HKL files to compare.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
	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
641
642
	if ( sym_str == NULL ) {
		sym_str = strdup("1");
643
	}
Thomas White's avatar
Thomas White committed
644
645
	sym = get_pointgroup(sym_str);
	free(sym_str);
646

647
648
649
	afile = strdup(argv[optind++]);
	bfile = strdup(argv[optind]);

650
	if ( pdb == NULL ) {
Thomas White's avatar
Thomas White committed
651
652
		ERROR("You must provide a PDB file with the unit cell.\n");
		exit(1);
653
654
655
656
657
	}

	cell = load_cell_from_pdb(pdb);
	free(pdb);

658
659
	list1_raw = read_reflections(afile);
	if ( list1_raw == NULL ) {
Thomas White's avatar
Thomas White committed
660
		ERROR("Couldn't read file '%s'\n", afile);
661
662
		return 1;
	}
Thomas White's avatar
Thomas White committed
663

664
665
	list2_raw = read_reflections(bfile);
	if ( list2_raw == NULL ) {
Thomas White's avatar
Thomas White committed
666
		ERROR("Couldn't read file '%s'\n", bfile);
667
668
669
		return 1;
	}

Thomas White's avatar
Thomas White committed
670
	/* Check that the intensities have the correct symmetry */
671
	if ( check_list_symmetry(list1_raw, sym) ) {
Thomas White's avatar
Thomas White committed
672
		ERROR("The first input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
673
		      " have symmetry %s\n", symmetry_name(sym));
Thomas White's avatar
Thomas White committed
674
675
		return 1;
	}
676
	if ( check_list_symmetry(list2_raw, sym) ) {
Thomas White's avatar
Thomas White committed
677
		ERROR("The second input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
678
		      " have symmetry %s\n", symmetry_name(sym));
Thomas White's avatar
Thomas White committed
679
680
681
		return 1;
	}

Thomas White's avatar
Thomas White committed
682
683
684
685
686
687
688
689
690
691
692
693
	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);

694
695
696
	list1 = asymmetric_indices(list1_raw, sym);
	list2 = asymmetric_indices(list2_raw, sym);

Thomas White's avatar
Thomas White committed
697
	/* Select reflections to be used */
698
	ncom = 0;
699
	nrej = 0;
Thomas White's avatar
Thomas White committed
700
701
702
703
	nneg = 0;
	nres = 0;
	list1_acc = reflist_new();
	list2_acc = reflist_new();
Thomas White's avatar
Thomas White committed
704
705
	for ( refl1 = first_refl(list1, &iter);
	      refl1 != NULL;
Thomas White's avatar
Thomas White committed
706
707
	      refl1 = next_refl(refl1, iter) )
	{
Thomas White's avatar
Thomas White committed
708
		signed int h, k, l;
709
		double val1, val2;
710
		double esd1, esd2;
Thomas White's avatar
Thomas White committed
711
		Reflection *refl2;
Thomas White's avatar
Thomas White committed
712
713
		Reflection *refl1_acc;
		Reflection *refl2_acc;
714

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

717
718
		refl2 = find_refl(list2, h, k, l);
		if ( refl2 == NULL ) continue;
719

Thomas White's avatar
Thomas White committed
720
721
722
		val1 = get_intensity(refl1);
		val2 = get_intensity(refl2);

723
724
725
726
727
728
729
730
731
732
		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
733
734
735
736
737
738
739
740
741
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
770
771
772
773
774
		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);
			if ( res > rmin_fix ) {
				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);

775
776
		ncom++;

777
	}
Thomas White's avatar
Thomas White committed
778

Thomas White's avatar
Thomas White committed
779
780
781
782
783
	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);
784
785
	}

Thomas White's avatar
Thomas White committed
786
787
788
789
	if ( config_ignorenegs && (nneg > 0) ) {
		STATUS("Discarded %i reflection pairs because either or both"
		       " versions had negative intensities.", nneg);
	}
790

Thomas White's avatar
Thomas White committed
791
792
793
	if ( config_zeronegs && (nneg > 0) ) {
		STATUS("For %i reflection pairs, either or both versions had"
		       " negative intensities which were set to zero.", nneg);
794
795
	}

Thomas White's avatar
Thomas White committed
796
797
798
	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
799
800
	}

Thomas White's avatar
Thomas White committed
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
	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,
	       nshells);

	reflist_free(list1_acc);
	reflist_free(list2_acc);

820
821
	return 0;
}