compare_hkl.c 21 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
}


318
struct shells
Thomas White's avatar
Thomas White committed
319
{
320
321
	int config_intshells;
	int nshells;
Thomas White's avatar
Thomas White committed
322
323
	double *rmins;
	double *rmaxs;
324
325
326
327
328
329
330
};


static struct shells *set_intensity_shells(double min_I, double max_I,
                                           int nshells)
{
	struct shells *s;
Thomas White's avatar
Thomas White committed
331
	int i;
Thomas White's avatar
Thomas White committed
332

333
334
335
336
	if ( min_I >= max_I ) {
		ERROR("Invalid intensity range.\n");
		return NULL;
	}
Thomas White's avatar
Thomas White committed
337

338
339
340
341
342
343
344
345
346
347
348
349
350
351
	/* Adjust minimum and maximum intensities to get the most densely
	 * populated part of the reflections */
	max_I = min_I + (max_I-min_I)/5000.0;

	s = malloc(sizeof(struct shells));
	if ( s == NULL ) return NULL;

	s->rmins = malloc(nshells*sizeof(double));
	s->rmaxs = malloc(nshells*sizeof(double));

	if ( (s->rmins==NULL) || (s->rmaxs==NULL) ) {
		ERROR("Couldn't allocate memory for shells.\n");
		free(s);
		return NULL;
Thomas White's avatar
Thomas White committed
352
353
	}

354
355
356
	s->config_intshells = 1;
	s->nshells = nshells;

Thomas White's avatar
Thomas White committed
357
	for ( i=0; i<nshells; i++ ) {
358
359
360
361

		s->rmins[i] = min_I + i*(max_I - min_I)/nshells;;
		s->rmaxs[i] = min_I + (i+1)*(max_I - min_I)/nshells;;

Thomas White's avatar
Thomas White committed
362
363
	}

364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
	return s;
}


static struct shells *set_resolution_shells(double rmin, double rmax,
                                            int nshells)
{
	struct shells *s;
	double total_vol, vol_per_shell;
	int i;

	s = malloc(sizeof(struct shells));
	if ( s == NULL ) return NULL;

	s->rmins = malloc(nshells*sizeof(double));
	s->rmaxs = malloc(nshells*sizeof(double));

	if ( (s->rmins==NULL) || (s->rmaxs==NULL) ) {
		ERROR("Couldn't allocate memory for resolution shells.\n");
		free(s);
		return NULL;
Thomas White's avatar
Thomas White committed
385
	}
386

387
388
389
	s->config_intshells = 0;
	s->nshells = nshells;

390
	total_vol = pow(rmax, 3.0) - pow(rmin, 3.0);
Thomas White's avatar
Thomas White committed
391
	vol_per_shell = total_vol / nshells;
392
	s->rmins[0] = rmin;
Thomas White's avatar
Thomas White committed
393
	for ( i=1; i<nshells; i++ ) {
394
395
396

		double r;

397
		r = vol_per_shell + pow(s->rmins[i-1], 3.0);
398
399
		r = pow(r, 1.0/3.0);

400
		/* Shells of constant volume */
401
402
		s->rmaxs[i-1] = r;
		s->rmins[i] = r;
403

404
405
406
407
		/* Shells of constant thickness */
		//rmins[i] = rmins[i-1] + (rmax-rmin)/NBINS;
		//rmaxs[i-1] = rmins[i-1] + (rmax-rmin)/NBINS;

408
	}
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
	s->rmaxs[nshells-1] = rmax;

	return s;
}


static double shell_label(struct shells *s, int i)
{
	if ( s->config_intshells ) {
		return (i+0.5) / s->nshells;
	} else {
		return s->rmins[i] + (s->rmaxs[i] - s->rmins[i])/2.0;
	}
}


static int get_bin(struct shells *s, Reflection *refl, UnitCell *cell)
{
	if ( s->config_intshells ) {

		double intensity;
		int bin, j;

		intensity = get_intensity(refl);

		bin = -1;
		for ( j=0; j<s->nshells; j++ ) {
			if ( (intensity>s->rmins[j])
			  && (intensity<=s->rmaxs[j]) )
			{
				bin = j;
				break;
			}
		}

		return bin;


	} else {
448

449
		double d;
450
451
		int bin, j;
		signed int h, k, l;
Thomas White's avatar
Thomas White committed
452

453
		get_indices(refl, &h, &k, &l);
454
		d = 2.0 * resolution(cell, h, k, l);
Thomas White's avatar
Thomas White committed
455

456
		bin = -1;
457
458
		for ( j=0; j<s->nshells; j++ ) {
			if ( (d>s->rmins[j]) && (d<=s->rmaxs[j]) ) {
459
460
461
462
				bin = j;
				break;
			}
		}
463

Thomas White's avatar
Thomas White committed
464
		/* Allow for slight rounding errors */
465
		if ( (bin == -1) && (d <= s->rmins[0]) ) bin = 0;
Thomas White's avatar
Thomas White committed
466
		assert(bin != -1);
Thomas White's avatar
Thomas White committed
467

468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
		return bin;

	}
}


static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
                   double rmin, double rmax, enum fom fom,
                   int config_unity, int nshells, const char *filename,
                   int config_intshells, double min_I, double max_I)
{
	int *cts;
	int i;
	Reflection *refl1;
	RefListIterator *iter;
	FILE *fh;
	double scale;
	struct fom_context *fctx;
	struct shells *shells;
	const char *t1, *t2;
	int n_out;

	cts = malloc(nshells*sizeof(int));
	fctx = init_fom(fom, num_reflections(list1), nshells);

	if ( (fctx==NULL) || (cts==NULL) ) {
		ERROR("Couldn't allocate memory for resolution shells.\n");
		return;
	}

	for ( i=0; i<nshells; i++ ) {
		cts[i] = 0;
	}

	if ( config_unity ) {
		scale = 1.0;
	} else {
		scale = stat_scale_intensity(list1, list2);
	}

	/* Calculate the bins */
	if ( config_intshells ) {
		shells = set_intensity_shells(min_I, max_I, nshells);
	} else {
		shells = set_resolution_shells(rmin, rmax, nshells);
	}

	if ( shells == NULL ) {
		ERROR("Failed to set up shells.\n");
		return;
	}

	n_out = 0;
	for ( refl1 = first_refl(list1, &iter);
	      refl1 != NULL;
	      refl1 = next_refl(refl1, iter) )
	{
		signed int h, k, l;
		int bin;
		double i1, i2;
		Reflection *refl2;

		get_indices(refl1, &h, &k, &l);

532
533
534
		refl2 = find_refl(list2, h, k, l);
		if ( refl2 == NULL ) continue;

535
536
537
538
539
540
		bin = get_bin(shells, refl1, cell);
		if ( bin == -1 ) {
			n_out++;
			continue;
		}

Thomas White's avatar
Thomas White committed
541
		i1 = get_intensity(refl1);
Thomas White's avatar
Thomas White committed
542
		i2 = scale * get_intensity(refl2);
543

Thomas White's avatar
Thomas White committed
544
545
546
		add_to_fom(fctx, i1, i2, bin);
		cts[bin]++;
	}
547
548
549
	if ( n_out)  {
		ERROR("Warning: %i reflection pairs outside range.\n", n_out);
	}
550

Thomas White's avatar
Thomas White committed
551
	switch ( fom ) {
552

Thomas White's avatar
Thomas White committed
553
554
555
		case FOM_R1I :
		STATUS("Overall R1(I) = %.2f %%\n", 100.0*fom_overall(fctx));
		break;
556

Thomas White's avatar
Thomas White committed
557
558
559
		case FOM_R1F :
		STATUS("Overall R1(F) = %.2f %%\n", 100.0*fom_overall(fctx));
		break;
560

Thomas White's avatar
Thomas White committed
561
562
563
		case FOM_R2 :
		STATUS("Overall R(2) = %.2f %%\n", 100.0*fom_overall(fctx));
		break;
564

Thomas White's avatar
Thomas White committed
565
566
567
		case FOM_RSPLIT :
		STATUS("Overall Rsplit = %.2f %%\n", 100.0*fom_overall(fctx));
		break;
Thomas White's avatar
Thomas White committed
568

Thomas White's avatar
Thomas White committed
569
570
571
572
573
574
575
		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
576

577
578
	}

579
	fh = fopen(filename, "w");
580
	if ( fh == NULL ) {
581
		ERROR("Couldn't open '%s'\n", filename);
582
583
584
		return;
	}

585
586
587
588
589
590
591
592
	if ( config_intshells ) {
		t1 = "Relative I  ";
		t2 = "";
	} else {
		t1 = "  1/d centre";
		t2 = "      d / A";
	}

Thomas White's avatar
Thomas White committed
593
594
595
	switch ( fom ) {

		case FOM_R1I :
596
		fprintf(fh, "%s  R1(I)/%%       nref%s\n", t1, t2);
Thomas White's avatar
Thomas White committed
597
		break;
598

Thomas White's avatar
Thomas White committed
599
		case FOM_R1F :
600
		fprintf(fh, "%s  R1(F)/%%       nref%s\n", t1, t2);
601
602
		break;

Thomas White's avatar
Thomas White committed
603
		case FOM_R2 :
604
		fprintf(fh, "%s     R2/%%       nref%s\n", t1, t2);
605
606
		break;

Thomas White's avatar
Thomas White committed
607
		case FOM_RSPLIT :
608
		fprintf(fh, "%s Rsplit/%%       nref%s\n", t1, t2);
609
610
		break;

Thomas White's avatar
Thomas White committed
611
		case FOM_CC :
612
		fprintf(fh, "%s       CC       nref%s\n", t1, t2);
Thomas White's avatar
Thomas White committed
613
614
615
		break;

		case FOM_CCSTAR :
616
		fprintf(fh, "%s      CC*       nref%s\n", t1, t2);
617
618
619
		break;

	}
620

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

623
		double r, cen;
624

625
		cen = shell_label(shells, i);
Thomas White's avatar
Thomas White committed
626
		r = fom_shell(fctx, i);
627

Thomas White's avatar
Thomas White committed
628
		switch ( fom ) {
629

Thomas White's avatar
Thomas White committed
630
631
632
633
		case FOM_R1I :
		case FOM_R1F :
		case FOM_R2 :
		case FOM_RSPLIT :
634
635
636
637
638
639
640
		if ( config_intshells ) {
			fprintf(fh, "%10.3f %10.2f %10i\n",
				cen, r*100.0, cts[i]);
		} else {
			fprintf(fh, "%10.3f %10.2f %10i %10.2f\n",
				cen*1.0e-9, r*100.0, cts[i], (1.0/cen)*1e10);
		}
Thomas White's avatar
Thomas White committed
641
		break;
642

Thomas White's avatar
Thomas White committed
643
644
		case FOM_CC :
		case FOM_CCSTAR :
645
646
647
648
649
650
651
		if ( config_intshells ) {
			fprintf(fh, "%10.3f %10.7f %10i\n",
				cen, r, cts[i]);
		} else {
			fprintf(fh, "%10.3f %10.7f %10i %10.2f\n",
				cen*1.0e-9, r, cts[i], (1.0/cen)*1e10);
		}
Thomas White's avatar
Thomas White committed
652
		break;
Thomas White's avatar
Thomas White committed
653
654
655

	}

Thomas White's avatar
Thomas White committed
656
	}
Thomas White's avatar
Thomas White committed
657

Thomas White's avatar
Thomas White committed
658
	fclose(fh);
Thomas White's avatar
Thomas White committed
659
660
661
}


662
663
664
int main(int argc, char *argv[])
{
	int c;
665
	UnitCell *cell;
666
667
	char *afile = NULL;
	char *bfile = NULL;
Thomas White's avatar
Thomas White committed
668
669
	char *sym_str = NULL;
	SymOpList *sym;
Thomas White's avatar
Thomas White committed
670
671
672
	int ncom, nrej, nneg, nres;
	RefList *list1_acc;
	RefList *list2_acc;
Thomas White's avatar
Thomas White committed
673
674
	RefList *list1;
	RefList *list2;
675
676
	RefList *list1_raw;
	RefList *list2_raw;
Thomas White's avatar
Thomas White committed
677
	enum fom fom = FOM_R1I;
678
	char *pdb = NULL;
679
680
	float rmin_fix = -1.0;
	float rmax_fix = -1.0;
Thomas White's avatar
Thomas White committed
681
	double rmin, rmax;
Thomas White's avatar
Thomas White committed
682
683
	Reflection *refl1;
	RefListIterator *iter;
684
	float sigma_cutoff = -INFINITY;
Thomas White's avatar
Thomas White committed
685
686
687
	int config_ignorenegs = 0;
	int config_zeronegs = 0;
	int config_unity = 0;
688
	int config_intshells = 0;
Thomas White's avatar
Thomas White committed
689
	int nshells = 10;
690
	char *shell_file = NULL;
691
692
	double min_I = +INFINITY;
	double max_I = -INFINITY;
693
694
695
696

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
697
		{"symmetry",           1, NULL,               'y'},
698
		{"pdb",                1, NULL,               'p'},
699
700
		{"rmin",               1, NULL,                2},
		{"rmax",               1, NULL,                3},
Thomas White's avatar
Thomas White committed
701
		{"fom",                1, NULL,                4},
702
		{"sigma-cutoff",       1, NULL,                5},
Thomas White's avatar
Thomas White committed
703
		{"nshells",            1, NULL,                6},
704
		{"shell-file",         1, NULL,                7},
Thomas White's avatar
Thomas White committed
705
706
		{"ignore-negs",        0, &config_ignorenegs,  1},
		{"zero-negs",          0, &config_zeronegs,    1},
707
		{"intensity-shells",   0, &config_intshells,   1},
708
709
710
711
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
712
	while ((c = getopt_long(argc, argv, "hy:p:u",
713
714
	                        longopts, NULL)) != -1)
	{
715
716

		switch (c) {
Thomas White's avatar
Thomas White committed
717
718

			case 'h' :
719
720
721
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
722
			case 'y' :
Thomas White's avatar
Thomas White committed
723
			sym_str = strdup(optarg);
724
725
			break;

Thomas White's avatar
Thomas White committed
726
			case 'p' :
727
728
729
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
730
			case 'u' :
731
732
733
			config_unity = 1;
			break;

Thomas White's avatar
Thomas White committed
734
			case 0 :
735
736
			break;

Thomas White's avatar
Thomas White committed
737
			case 2 :
738
739
740
741
742
743
			if ( sscanf(optarg, "%e", &rmin_fix) != 1 ) {
				ERROR("Invalid value for --rmin\n");
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
744
			case 3 :
745
746
747
748
749
750
			if ( sscanf(optarg, "%e", &rmax_fix) != 1 ) {
				ERROR("Invalid value for --rmax\n");
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
751
			case 4 :
Thomas White's avatar
Thomas White committed
752
			fom = get_fom(optarg);
753
754
			break;

755
756
757
758
759
760
761
			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
762
763
764
765
766
767
768
			case 6 :
			if ( sscanf(optarg, "%i", &nshells) != 1 ) {
				ERROR("Invalid value for --nshells\n");
				return 1;
			}
			break;

769
770
771
772
			case 7 :
			shell_file = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
773
			default :
774
775
			ERROR("Unhandled option '%c'\n", c);
			break;
Thomas White's avatar
Thomas White committed
776

777
778
779
780
		}

	}

781
782
783
784
785
	if ( argc != (optind+2) ) {
		ERROR("Please provide exactly two HKL files to compare.\n");
		return 1;
	}

Thomas White's avatar
Thomas White committed
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
	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
807
808
	if ( sym_str == NULL ) {
		sym_str = strdup("1");
809
	}
Thomas White's avatar
Thomas White committed
810
811
	sym = get_pointgroup(sym_str);
	free(sym_str);
812

813
814
815
	afile = strdup(argv[optind++]);
	bfile = strdup(argv[optind]);

816
	if ( pdb == NULL ) {
Thomas White's avatar
Thomas White committed
817
818
		ERROR("You must provide a PDB file with the unit cell.\n");
		exit(1);
819
820
	}

821
822
	if ( shell_file == NULL ) shell_file = strdup("shells.dat");

823
824
825
	cell = load_cell_from_pdb(pdb);
	free(pdb);

826
827
	list1_raw = read_reflections(afile);
	if ( list1_raw == NULL ) {
Thomas White's avatar
Thomas White committed
828
		ERROR("Couldn't read file '%s'\n", afile);
829
830
		return 1;
	}
Thomas White's avatar
Thomas White committed
831

832
833
	list2_raw = read_reflections(bfile);
	if ( list2_raw == NULL ) {
Thomas White's avatar
Thomas White committed
834
		ERROR("Couldn't read file '%s'\n", bfile);
835
836
837
		return 1;
	}

Thomas White's avatar
Thomas White committed
838
	/* Check that the intensities have the correct symmetry */
839
	if ( check_list_symmetry(list1_raw, sym) ) {
Thomas White's avatar
Thomas White committed
840
		ERROR("The first input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
841
		      " have symmetry %s\n", symmetry_name(sym));
Thomas White's avatar
Thomas White committed
842
843
		return 1;
	}
844
	if ( check_list_symmetry(list2_raw, sym) ) {
Thomas White's avatar
Thomas White committed
845
		ERROR("The second input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
846
		      " have symmetry %s\n", symmetry_name(sym));
Thomas White's avatar
Thomas White committed
847
848
849
		return 1;
	}

Thomas White's avatar
Thomas White committed
850
851
852
853
854
855
856
857
858
859
860
861
	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);

862
863
864
	list1 = asymmetric_indices(list1_raw, sym);
	list2 = asymmetric_indices(list2_raw, sym);

Thomas White's avatar
Thomas White committed
865
	/* Select reflections to be used */
866
	ncom = 0;
867
	nrej = 0;
Thomas White's avatar
Thomas White committed
868
869
870
871
	nneg = 0;
	nres = 0;
	list1_acc = reflist_new();
	list2_acc = reflist_new();
Thomas White's avatar
Thomas White committed
872
873
	for ( refl1 = first_refl(list1, &iter);
	      refl1 != NULL;
Thomas White's avatar
Thomas White committed
874
875
	      refl1 = next_refl(refl1, iter) )
	{
Thomas White's avatar
Thomas White committed
876
		signed int h, k, l;
877
		double val1, val2;
878
		double esd1, esd2;
Thomas White's avatar
Thomas White committed
879
		Reflection *refl2;
Thomas White's avatar
Thomas White committed
880
881
		Reflection *refl1_acc;
		Reflection *refl2_acc;
882

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

885
886
		refl2 = find_refl(list2, h, k, l);
		if ( refl2 == NULL ) continue;
887

Thomas White's avatar
Thomas White committed
888
889
890
		val1 = get_intensity(refl1);
		val2 = get_intensity(refl2);

891
892
893
894
895
896
897
898
899
900
		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
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
		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
929
			if ( res > rmax_fix ) {
Thomas White's avatar
Thomas White committed
930
931
932
933
934
935
936
937
938
939
940
941
942
				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);

943
944
945
		if ( val1 > max_I ) max_I = val1;
		if ( val1 < min_I ) min_I = val1;

946
947
		ncom++;

948
	}
Thomas White's avatar
Thomas White committed
949

Thomas White's avatar
Thomas White committed
950
951
952
953
954
	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);
955
956
	}

Thomas White's avatar
Thomas White committed
957
958
	if ( config_ignorenegs && (nneg > 0) ) {
		STATUS("Discarded %i reflection pairs because either or both"
Thomas White's avatar
Thomas White committed
959
		       " versions had negative intensities.\n", nneg);
Thomas White's avatar
Thomas White committed
960
	}
961

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

Thomas White's avatar
Thomas White committed
967
968
969
	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
970
971
	}

Thomas White's avatar
Thomas White committed
972
973
974
975
976
977
978
979
980
981
982
983
984
	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);

985
986
987
988
989
990
991
992
993
994
995
	if ( rmin_fix >= 0.0 ) {
		rmin = rmin_fix;
	}
	if ( rmax_fix >= 0.0 ) {
		rmax = rmax_fix;
	}
	if ( (rmin_fix>=0.0) || (rmax_fix>=0.0) ) {
		STATUS("Fixed resolution range: %f to %f nm^-1"
		       " (%.2f to %.2f Angstroms).\n",
		       rmin/1e9, rmax/1e9, 1e10/rmin, 1e10/rmax);
	}
Thomas White's avatar
Thomas White committed
996
	do_fom(list1_acc, list2_acc, cell, rmin, rmax, fom, config_unity,
997
	       nshells, shell_file, config_intshells, min_I, max_I);
Thomas White's avatar
Thomas White committed
998

999
	free(shell_file);
Thomas White's avatar
Thomas White committed
1000
	reflist_free(list1_acc);