aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-03-16 18:23:48 +0100
committerThomas White <taw@physics.org>2012-02-22 15:27:20 +0100
commit86e621ad289638d7bf1440092cffd4c2ac64f982 (patch)
tree0abc051952d4d7d581bfe7b2d7403972ebd7327f /src
parenteaee0f6415b06f40ff7f477328f9f26f01c4d7ec (diff)
Remove old reflections.h
Diffstat (limited to 'src')
-rw-r--r--src/check_hkl.c92
-rw-r--r--src/compare_hkl.c220
-rw-r--r--src/get_hkl.c283
-rw-r--r--src/indexamajig.c2
-rw-r--r--src/partialator.c4
-rw-r--r--src/pattern_sim.c30
-rw-r--r--src/povray.c47
-rw-r--r--src/povray.h8
-rw-r--r--src/process_hkl.c9
-rw-r--r--src/reflections.c230
-rw-r--r--src/reflections.h38
-rw-r--r--src/reflist-utils.c247
-rw-r--r--src/reflist-utils.h18
-rw-r--r--src/reflist.c27
-rw-r--r--src/reflist.h6
-rw-r--r--src/render_hkl.c53
-rw-r--r--src/statistics.c385
-rw-r--r--src/statistics.h48
-rw-r--r--src/stream.c49
19 files changed, 927 insertions, 869 deletions
diff --git a/src/check_hkl.c b/src/check_hkl.c
index fb6d02d5..903188a1 100644
--- a/src/check_hkl.c
+++ b/src/check_hkl.c
@@ -22,9 +22,10 @@
#include <getopt.h>
#include "utils.h"
-#include "reflections.h"
#include "statistics.h"
#include "symmetry.h"
+#include "reflist.h"
+#include "reflist-utils.h"
/* Number of bins for plot of resolution shells */
@@ -46,9 +47,8 @@ static void show_help(const char *s)
}
-static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell,
- const char *sym, unsigned int *counts,
- const double *sigma, double rmin_fix, double rmax_fix)
+static void plot_shells(RefList *list, UnitCell *cell, const char *sym,
+ double rmin_fix, double rmax_fix)
{
double num[NBINS];
int cts[NBINS];
@@ -70,6 +70,8 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell,
int nmeas = 0;
int nmeastot = 0;
int nout = 0;
+ Reflection *refl;
+ RefListIterator *iter;
if ( cell == NULL ) {
ERROR("Need the unit cell to plot resolution shells.\n");
@@ -93,16 +95,17 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell,
mean[i] = 0;
}
- rmin = +INFINITY;
- rmax = 0.0;
- for ( i=0; i<num_items(items); i++ ) {
+ /* Iterate over all common reflections and calculate min and max
+ * resolution */
+ rmin = +INFINITY; rmax = 0.0;
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
- struct refl_item *it;
signed int h, k, l;
double d;
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
+ get_indices(refl, &h, &k, &l);
d = resolution(cell, h, k, l) * 2.0;
if ( d > rmax ) rmax = d;
@@ -179,16 +182,16 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell,
free(counted);
/* Calculate means */
- for ( i=0; i<num_items(items); i++ ) {
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
- struct refl_item *it;
signed int h, k, l;
double d;
int bin;
int j;
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
+ get_indices(refl, &h, &k, &l);
d = resolution(cell, h, k, l) * 2.0;
@@ -205,7 +208,7 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell,
}
measured[bin]++;
- mean[bin] += lookup_intensity(ref, h, k, l);
+ mean[bin] += get_intensity(refl);
}
@@ -214,21 +217,21 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell,
}
/* Characterise the data set */
- for ( i=0; i<num_items(items); i++ ) {
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
- struct refl_item *it;
signed int h, k, l;
double d;
int bin;
int j;
double val, esd;
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
+ get_indices(refl, &h, &k, &l);
d = resolution(cell, h, k, l) * 2.0;
- val = lookup_intensity(ref, h, k, l);
- esd = lookup_sigma(sigma, h, k, l);
+ val = get_intensity(refl);
+ esd = get_esd_intensity(refl);
bin = -1;
for ( j=0; j<NBINS; j++ ) {
@@ -245,11 +248,11 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell,
if ( !isfinite(val/esd) ) continue;
/* measured[bin] was done earlier */
- measurements[bin] += lookup_count(counts, h, k, l);
+ measurements[bin] += get_redundancy(refl);
snr[bin] += val / esd;
snr_total += val / esd;
nmeas++;
- nmeastot += lookup_count(counts, h, k, l);
+ nmeastot += get_redundancy(refl);
var[bin] += pow(val-mean[bin], 2.0);
@@ -289,17 +292,15 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell,
int main(int argc, char *argv[])
{
int c;
- double *ref;
UnitCell *cell;
char *file = NULL;
char *sym = NULL;
- int i;
- ReflItemList *items;
- ReflItemList *good_items;
+ RefList *raw_list;
+ RefList *list;
+ Reflection *refl;
+ RefListIterator *iter;
char *pdb = NULL;
- double *esd;
int rej = 0;
- unsigned int *cts;
float rmin_fix = -1.0;
float rmax_fix = -1.0;
@@ -369,37 +370,35 @@ int main(int argc, char *argv[])
cell = load_cell_from_pdb(pdb);
free(pdb);
- ref = new_list_intensity();
- esd = new_list_sigma();
- cts = new_list_count();
- items = read_reflections(file, ref, NULL, cts, esd);
- if ( items == NULL ) {
- ERROR("Couldn't open file '%s'\n", file);
+ raw_list = read_reflections(file);
+ if ( raw_list == NULL ) {
+ ERROR("Couldn't read file '%s'\n", file);
return 1;
}
/* Check that the intensities have the correct symmetry */
- if ( check_symmetry(items, sym) ) {
+ if ( check_list_symmetry(raw_list, sym) ) {
ERROR("The input reflection list does not appear to"
" have symmetry %s\n", sym);
return 1;
}
- /* Reject reflections */
- good_items = new_items();
- for ( i=0; i<num_items(items); i++ ) {
+ /* Reject some reflections */
+ list = reflist_new();
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
- struct refl_item *it;
signed int h, k, l;
double val, sig;
int ig = 0;
double d;
+ Reflection *new;
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
+ get_indices(refl, &h, &k, &l);
- val = lookup_intensity(ref, h, k, l);
- sig = lookup_sigma(esd, h, k, l);
+ val = get_intensity(refl);
+ sig = get_esd_intensity(refl);
if ( val < 3.0 * sig ) {
rej++;
@@ -412,11 +411,12 @@ int main(int argc, char *argv[])
//if ( ig ) continue;
- add_item(good_items, h, k, l);
+ new = add_refl(list, h, k, l);
+ copy_data(new, refl);
}
- plot_shells(ref, items, cell, sym, cts, esd, rmin_fix, rmax_fix);
+ plot_shells(list, cell, sym, rmin_fix, rmax_fix);
return 0;
}
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index abe8bc6a..16a744ea 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -20,11 +20,12 @@
#include <string.h>
#include <unistd.h>
#include <getopt.h>
+#include <assert.h>
#include "utils.h"
-#include "reflections.h"
#include "statistics.h"
#include "symmetry.h"
+#include "reflist-utils.h"
/* Number of bins for plot of resolution shells */
@@ -38,7 +39,7 @@ static void show_help(const char *s)
"Compare intensity lists.\n"
"\n"
" -h, --help Display this help message.\n"
-" -o, --output=<filename> Specify output filename for correction factor.\n"
+" -o, --ratio=<filename> Specify output filename for ratios.\n"
" -y, --symmetry=<sym> The symmetry of both the input files.\n"
" -p, --pdb=<filename> PDB file to use (default: molecule.pdb).\n"
" --shells Plot the figures of merit by resolution.\n"
@@ -48,9 +49,9 @@ static void show_help(const char *s)
}
-static void plot_shells(const double *ref1, const double *ref2,
- ReflItemList *items, double scale, UnitCell *cell,
- const char *sym, double rmin_fix, double rmax_fix)
+static void plot_shells(RefList *list1, RefList *list2, double scale,
+ UnitCell *cell, const char *sym,
+ double rmin_fix, double rmax_fix)
{
double num[NBINS];
int cts[NBINS];
@@ -69,6 +70,12 @@ static void plot_shells(const double *ref1, const double *ref2,
int ctot;
FILE *fh;
int nout = 0;
+ Reflection *refl1;
+ RefListIterator *iter;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ int hmax, kmax, lmax;
if ( cell == NULL ) {
ERROR("Need the unit cell to plot resolution shells.\n");
@@ -90,16 +97,20 @@ static void plot_shells(const double *ref1, const double *ref2,
snr[i] = 0;
}
- rmin = +INFINITY;
- rmax = 0.0;
- for ( i=0; i<num_items(items); i++ ) {
+ /* Iterate over all common reflections and calculate min and max
+ * resolution */
+ rmin = +INFINITY; rmax = 0.0;
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) ) {
- struct refl_item *it;
signed int h, k, l;
double d;
+ Reflection *refl2;
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
d = resolution(cell, h, k, l) * 2.0;
if ( d > rmax ) rmax = d;
@@ -117,6 +128,7 @@ static void plot_shells(const double *ref1, const double *ref2,
if ( rmin_fix > 0.0 ) rmin = rmin_fix;
if ( rmax_fix > 0.0 ) rmax = rmax_fix;
+ /* Calculate the resolution bins */
total_vol = pow(rmax, 3.0) - pow(rmin, 3.0);
vol_per_shell = total_vol / NBINS;
rmins[0] = rmin;
@@ -145,9 +157,16 @@ static void plot_shells(const double *ref1, const double *ref2,
/* Count the number of reflections possible in each shell */
counted = new_list_count();
- for ( h=-50; h<=+50; h++ ) {
- for ( k=-50; k<=+50; k++ ) {
- for ( l=-50; l<=+50; l++ ) {
+ cell_get_reciprocal(cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+ hmax = rmax / modulus(asx, asy, asz);
+ kmax = rmax / modulus(bsx, bsy, bsz);
+ lmax = rmax / modulus(csx, csy, csz);
+
+ for ( h=-hmax; h<hmax; h++ ) {
+ for ( k=-kmax; k<kmax; k++ ) {
+ for ( l=-lmax; l<lmax; l++ ) {
double d;
signed int hs, ks, ls;
@@ -178,17 +197,21 @@ static void plot_shells(const double *ref1, const double *ref2,
den = 0.0;
ctot = 0;
nout = 0;
- for ( i=0; i<num_items(items); i++ ) {
- struct refl_item *it;
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) ) {
+
signed int h, k, l;
double d;
int bin;
double i1, i2;
int j;
+ Reflection *refl2;
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
d = resolution(cell, h, k, l) * 2.0;
@@ -206,8 +229,8 @@ static void plot_shells(const double *ref1, const double *ref2,
continue;
}
- i1 = lookup_intensity(ref1, h, k, l);
- i2 = lookup_intensity(ref2, h, k, l);
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
i2 *= scale;
num[bin] += fabs(i1 - i2);
@@ -238,33 +261,32 @@ static void plot_shells(const double *ref1, const double *ref2,
int main(int argc, char *argv[])
{
int c;
- double *ref1;
- double *ref2;
- double *ref2_transformed;
- double *out;
UnitCell *cell;
- char *outfile = NULL;
+ char *ratiofile = NULL;
char *afile = NULL;
char *bfile = NULL;
char *sym = NULL;
double scale, scale_r2, scale_rdig, R1, R2, R1i, Rdiff, pearson;
double scale_rintint, scale_r1i, scale_r1, scale_r1fi;
- int i, ncom;
- ReflItemList *i1, *i2, *icommon;
+ int ncom;
+ RefList *list1;
+ RefList *list2;
+ RefList *list2_transformed;
+ RefList *ratio;
+ RefList *deleteme;
int config_shells = 0;
char *pdb = NULL;
- double *esd1;
- double *esd2;
int rej1 = 0;
int rej2 = 0;
- unsigned int *cts1;
float rmin_fix = -1.0;
float rmax_fix = -1.0;
+ Reflection *refl1;
+ RefListIterator *iter;
/* Long options */
const struct option longopts[] = {
{"help", 0, NULL, 'h'},
- {"output", 1, NULL, 'o'},
+ {"ratio" , 1, NULL, 'o'},
{"symmetry", 1, NULL, 'y'},
{"shells", 0, &config_shells, 1},
{"pdb", 1, NULL, 'p'},
@@ -282,7 +304,7 @@ int main(int argc, char *argv[])
return 0;
case 'o' :
- outfile = strdup(optarg);
+ ratiofile = strdup(optarg);
break;
case 'y' :
@@ -335,63 +357,62 @@ int main(int argc, char *argv[])
cell = load_cell_from_pdb(pdb);
free(pdb);
- ref1 = new_list_intensity();
- esd1 = new_list_sigma();
- cts1 = new_list_count();
- i1 = read_reflections(afile, ref1, NULL, cts1, esd1);
- if ( i1 == NULL ) {
- ERROR("Couldn't open file '%s'\n", afile);
+ list1 = read_reflections(afile);
+ if ( list1 == NULL ) {
+ ERROR("Couldn't read file '%s'\n", afile);
return 1;
}
- ref2 = new_list_intensity();
- esd2 = new_list_sigma();
- i2 = read_reflections(bfile, ref2, NULL, NULL, esd2);
- if ( i2 == NULL ) {
- ERROR("Couldn't open file '%s'\n", bfile);
+
+ list2 = read_reflections(bfile);
+ if ( list2 == NULL ) {
+ ERROR("Couldn't read file '%s'\n", bfile);
return 1;
}
/* Check that the intensities have the correct symmetry */
- if ( check_symmetry(i1, sym) ) {
+ if ( check_list_symmetry(list1, sym) ) {
ERROR("The first input reflection list does not appear to"
" have symmetry %s\n", sym);
return 1;
}
- if ( check_symmetry(i2, sym) ) {
+ if ( check_list_symmetry(list2, sym) ) {
ERROR("The second input reflection list does not appear to"
" have symmetry %s\n", sym);
return 1;
}
- /* List for output scale factor map */
- out = new_list_intensity();
-
/* Find common reflections (taking symmetry into account) */
- icommon = new_items();
- ref2_transformed = new_list_intensity();
- for ( i=0; i<num_items(i1); i++ ) {
+ list2_transformed = reflist_new();
+ ratio = reflist_new();
+ deleteme = reflist_new();
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) ) {
- struct refl_item *it;
signed int h, k, l;
signed int he, ke, le;
double val1, val2;
double sig1, sig2;
int ig = 0;
double d;
+ Reflection *refl2;
+ Reflection *tr;
- it = get_item(i1, i);
- h = it->h; k = it->k; l = it->l;
+ get_indices(refl1, &h, &k, &l);
- if ( !find_unique_equiv(i2, h, k, l, sym, &he, &ke, &le) ) {
- //STATUS("%i %i %i not matched (%f nm).\n", h, k, l,
- // 1.0/(2.0*resolution(cell, h, k, l)/1e9));
+ if ( !find_equiv_in_list(list2, h, k, l, sym, &he, &ke, &le) ) {
+ /* No common reflection */
+ add_refl(deleteme, h, k, l);
continue;
}
- val1 = lookup_intensity(ref1, h, k, l);
- val2 = lookup_intensity(ref2, he, ke, le);
- sig1 = lookup_sigma(esd1, h, k, l);
- sig2 = lookup_sigma(esd2, he, ke, le);
+ refl2 = find_refl(list2, he, ke, le);
+ assert(refl2 != NULL);
+
+ val1 = get_intensity(refl1);
+ val2 = get_intensity(refl2);
+ sig1 = get_esd_intensity(refl1);
+ sig2 = get_esd_intensity(refl2);
if ( val1 < 3.0 * sig1 ) {
rej1++;
@@ -408,67 +429,90 @@ int main(int argc, char *argv[])
//if ( ig ) continue;
- set_intensity(ref2_transformed, h, k, l, val2);
- set_intensity(out, h, k, l, val1/val2);
- add_item(icommon, h, k, l);
+ /* Add the old data from 'refl2' to a new list with the same
+ * indices as its equivalent in 'list1' */
+ tr = add_refl(list2_transformed, h, k, l);
+ copy_data(tr, refl2);
+
+ /* Add divided version to 'output' list */
+ tr = add_refl(ratio, h, k, l);
+ set_int(tr, val1/val2);
}
- ncom = num_items(icommon);
+
STATUS("%i reflections in '%s' had I < 3.0*sigma(I)\n", rej1, afile);
STATUS("%i reflections in '%s' had I < 3.0*sigma(I)\n", rej2, bfile);
+ ncom = num_reflections(list2_transformed);
STATUS("%i,%i reflections: %i in common\n",
- num_items(i1), num_items(i2), ncom);
+ num_reflections(list1), num_reflections(list2), ncom);
+
+ /* Trim reflections from 'list1' which had no equivalents in 'list2' */
+ for ( refl1 = first_refl(deleteme, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) ) {
- R1 = stat_r1_ignore(ref1, ref2_transformed, icommon, &scale_r1fi);
- STATUS("R1(F) = %5.4f %% (scale=%5.2e) (ignoring negative intensities)\n",
+ signed int h, k, l;
+ Reflection *del;
+
+ get_indices(refl1, &h, &k, &l);
+ del = find_refl(list1, h, k, l);
+ assert(del != NULL);
+
+ delete_refl(del);
+
+ }
+ reflist_free(deleteme);
+ reflist_free(list2);
+
+ R1 = stat_r1_ignore(list1, list2_transformed, &scale_r1fi);
+ STATUS("R1(F) = %5.4f %% (scale=%5.2e)"
+ " (ignoring negative intensities)\n",
R1*100.0, scale_r1fi);
- R1 = stat_r1_zero(ref1, ref2_transformed, icommon, &scale_r1);
- STATUS("R1(F) = %5.4f %% (scale=%5.2e) (zeroing negative intensities)\n",
+ R1 = stat_r1_zero(list1, list2_transformed, &scale_r1);
+ STATUS("R1(F) = %5.4f %% (scale=%5.2e)"
+ " (zeroing negative intensities)\n",
R1*100.0, scale_r1);
- R2 = stat_r2(ref1, ref2_transformed, icommon, &scale_r2);
+ R2 = stat_r2(list1, list2_transformed, &scale_r2);
STATUS("R2(I) = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale_r2);
- R1i = stat_r1_i(ref1, ref2_transformed, icommon, &scale_r1i);
+ R1i = stat_r1_i(list1, list2_transformed, &scale_r1i);
STATUS("R1(I) = %5.4f %% (scale=%5.2e)\n", R1i*100.0, scale_r1i);
- Rdiff = stat_rdiff_ignore(ref1, ref2_transformed, icommon, &scale_rdig);
- STATUS("Rint(F) = %5.4f %% (scale=%5.2e) (ignoring negative intensities)\n",
+ Rdiff = stat_rdiff_ignore(list1, list2_transformed, &scale_rdig);
+ STATUS("Rint(F) = %5.4f %% (scale=%5.2e)"
+ " (ignoring negative intensities)\n",
Rdiff*100.0, scale_rdig);
- Rdiff = stat_rdiff_zero(ref1, ref2_transformed, icommon, &scale);
- STATUS("Rint(F) = %5.4f %% (scale=%5.2e) (zeroing negative intensities)\n",
+ Rdiff = stat_rdiff_zero(list1, list2_transformed, &scale);
+ STATUS("Rint(F) = %5.4f %% (scale=%5.2e)"
+ " (zeroing negative intensities)\n",
Rdiff*100.0, scale);
- Rdiff = stat_rdiff_intensity(ref1, ref2_transformed, icommon,
- &scale_rintint);
+ Rdiff = stat_rdiff_intensity(list1, list2_transformed, &scale_rintint);
STATUS("Rint(I) = %5.4f %% (scale=%5.2e)\n",
Rdiff*100.0, scale_rintint);
- pearson = stat_pearson_i(ref1, ref2_transformed, icommon);
+ pearson = stat_pearson_i(list1, list2_transformed);
STATUS("Pearson r(I) = %5.4f\n", pearson);
- pearson = stat_pearson_f_ignore(ref1, ref2_transformed, icommon);
+ pearson = stat_pearson_f_ignore(list1, list2_transformed);
STATUS("Pearson r(F) = %5.4f (ignoring negative intensities)\n",
pearson);
- pearson = stat_pearson_f_zero(ref1, ref2_transformed, icommon);
+ pearson = stat_pearson_f_zero(list1, list2_transformed);
STATUS("Pearson r(F) = %5.4f (zeroing negative intensities)\n",
pearson);
if ( config_shells ) {
- plot_shells(ref1, ref2_transformed, icommon, scale_r1fi,
+ plot_shells(list1, list2_transformed, scale_r1fi,
cell, sym, rmin_fix, rmax_fix);
}
- if ( outfile != NULL ) {
-
- write_reflections(outfile, icommon, out, NULL,
- NULL, NULL, cell);
- STATUS("Sigma(I) values in output file are not meaningful.\n");
-
+ if ( ratiofile != NULL ) {
+ write_reflist(ratiofile, ratio, cell);
}
return 0;
diff --git a/src/get_hkl.c b/src/get_hkl.c
index 1b7e4406..026198ab 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -22,7 +22,7 @@
#include <getopt.h>
#include "utils.h"
-#include "reflections.h"
+#include "reflist-utils.h"
#include "symmetry.h"
#include "beam-parameters.h"
@@ -45,7 +45,8 @@ static void show_help(const char *s)
"To calculate Poisson samples accurately, you must also give:\n"
" -b, --beam=<file> Get beam parameters from file.\n"
"\n"
-"You can artificially 'twin' the reflections, or expand them out:\n"
+"You can artificially 'twin' the reflections, or expand them out. You can also"
+" do both, in which case the 'twinning' will be done first:\n"
" -w, --twin=<sym> Generate twinned data according to the given\n"
" point group.\n"
" -e, --expand=<sym> Expand reflections to this point group.\n"
@@ -59,99 +60,123 @@ static void show_help(const char *s)
"\n"
"Don't forget to specify the output filename:\n"
" -o, --output=<filename> Output filename (default: stdout).\n"
+" -p, --pdb=<filename> Use unit cell parameters from this PDB file to\n"
+" generate resolution values in the output file.\n"
);
}
/* Apply Poisson noise to all reflections */
-static void poisson_reflections(double *ref, ReflItemList *items,
- double adu_per_photon)
+static void poisson_reflections(RefList *list, double adu_per_photon)
{
- int i;
- const int n = num_items(items);
+ Reflection *refl;
+ RefListIterator *iter;
- for ( i=0; i<n; i++ ) {
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
- struct refl_item *it;
- double val;
- int c;
+ double val, c;
- it = get_item(items, i);
+ val = get_intensity(refl);
- val = lookup_intensity(ref, it->h, it->k, it->l);
c = adu_per_photon * poisson_noise(val/adu_per_photon);
- set_intensity(ref, it->h, it->k, it->l, c);
-
- progress_bar(i, n-1, "Simulating noise");
+ set_int(refl, c);
}
}
/* Apply 10% uniform noise to all reflections */
-static void noise_reflections(double *ref, ReflItemList *items)
+static void noise_reflections(RefList *list)
{
- int i;
- const int n = num_items(items);
-
- for ( i=0; i<n; i++ ) {
+ Reflection *refl;
+ RefListIterator *iter;
- struct refl_item *it;
- double val;
- double r;
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
- it = get_item(items, i);
+ double val, r;
- val = lookup_intensity(ref, it->h, it->k, it->l);
+ val = get_intensity(refl);
r = (double)random()/RAND_MAX;
val += 0.1 * val * r;
- set_intensity(ref, it->h, it->k, it->l, val);
-
- progress_bar(i, n-1, "Simulating noise");
+ set_int(refl, val);
}
}
-static ReflItemList *twin_reflections(double *ref, ReflItemList *items,
- const char *holo, const char *mero,
- double *esds)
+static RefList *template_reflections(RefList *list, RefList *template)
{
- int i;
- ReflItemList *new;
+ Reflection *refl;
+ RefListIterator *iter;
+ RefList *out;
- new = new_items();
+ out = reflist_new();
+
+ for ( refl = first_refl(template, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
+
+ signed int h, k, l;
+ Reflection *new;
+ Reflection *old;
+
+ get_indices(refl, &h, &k, &l);
+
+ old = find_refl(list, h, k, l);
+ if ( old == NULL ) continue;
+
+ new = add_refl(out, h, k, l);
+ copy_data(new, old);
+
+ }
+
+ return out;
+}
+
+
+static RefList *twin_reflections(RefList *in,
+ const char *holo, const char *mero)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+ RefList *out;
if ( num_general_equivs(holo) < num_general_equivs(mero) ) {
ERROR("%s is not a subgroup of %s!\n", mero, holo);
return NULL;
}
- for ( i=0; i<num_items(items); i++ ) {
+ out = reflist_new();
+
+ for ( refl = first_refl(in, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
double total, sigma;
- struct refl_item *it;
signed int h, k, l;
int n, j;
int skip;
- it = get_item(items, i);
+ get_indices(refl, &h, &k, &l);
/* There is a many-to-one correspondence between reflections
* in the merohedral and holohedral groups. Do the calculation
* only once for each reflection in the holohedral group, which
* contains fewer reflections.
*/
- get_asymm(it->h, it->k, it->l, &h, &k, &l, holo);
- if ( find_item(new, h, k, l) ) continue;
-
- n = num_equivs(h, k, l, holo);
+ get_asymm(h, k, l, &h, &k, &l, holo);
+ if ( find_refl(out, h, k, l) != NULL ) continue;
total = 0.0;
sigma = 0.0;
skip = 0;
+ n = num_equivs(h, k, l, holo);
for ( j=0; j<n; j++ ) {
signed int he, ke, le;
@@ -164,8 +189,8 @@ static ReflItemList *twin_reflections(double *ref, ReflItemList *items,
* equivalent which belongs to our definition of the
* asymmetric unit cell, so check them all.
*/
- if ( !find_unique_equiv(items, he, ke, le, mero,
- &hu, &ku, &lu) ) {
+ if ( !find_equiv_in_list(in, he, ke, le, mero,
+ &hu, &ku, &lu) ) {
/* Don't have this reflection, so bail out */
ERROR("Twinning %i %i %i requires the %i %i %i "
"reflection (or an equivalent in %s), "
@@ -176,61 +201,57 @@ static ReflItemList *twin_reflections(double *ref, ReflItemList *items,
break;
}
- total += lookup_intensity(ref, hu, ku, lu);
- sigma += pow(lookup_sigma(esds, hu, ku, lu), 2.0);
+ total += get_intensity(refl);
+ sigma += pow(get_esd_intensity(refl), 2.0);
}
if ( !skip ) {
- set_intensity(ref, h, k, l, total);
- set_sigma(esds, h, k, l, sqrt(sigma));
- add_item(new, h, k, l);
+ Reflection *new = add_refl(out, h, k, l);
+ set_int(new, total);
+ set_esd_intensity(new, sqrt(sigma));
}
}
- return new;
+ return out;
}
-static ReflItemList *expand_reflections(double *ref, ReflItemList *items,
- const char *target, const char *initial)
+static RefList *expand_reflections(RefList *in,
+ const char *target, const char *initial)
{
- int i;
- ReflItemList *new;
-
- new = new_items();
+ Reflection *refl;
+ RefListIterator *iter;
+ RefList *out;
if ( num_general_equivs(target) > num_general_equivs(initial) ) {
ERROR("%s is not a subgroup of %s!\n", initial, target);
return NULL;
}
- for ( i=0; i<num_items(items); i++ ) {
+ out = reflist_new();
+
+ for ( refl = first_refl(in, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
- struct refl_item *it;
signed int h, k, l;
- signed int hd, kd, ld;
int n, j;
double intensity;
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
-
- /* Actually we don't really care what the equivalent is,
- * we just want to be sure that there is nly be one version of
- * this reflection. */
- find_unique_equiv(items, h, k, l, initial, &hd, &kd, &ld);
+ get_indices(refl, &h, &k, &l);
+ intensity = get_intensity(refl);
- /* Now find out how many reflections need to be filled in */
n = num_equivs(h, k, l, initial);
- intensity = lookup_intensity(ref, h, k, l);
+ /* For each equivalent in the higher symmetry group */
for ( j=0; j<n; j++ ) {
signed int he, ke, le;
+ Reflection *new;
/* Get the equivalent */
get_equiv(h, k, l, &he, &ke, &le, initial, j);
@@ -239,43 +260,33 @@ static ReflItemList *expand_reflections(double *ref, ReflItemList *items,
get_asymm(he, ke, le, &he, &ke, &le, target);
/* Make sure the intensity is in the right place */
- set_intensity(ref, he, ke, le, intensity);
-
- /* Add the reflection, but only once */
- if ( !find_item(new, he, ke, le) ) {
- add_item(new, he, ke, le);
- }
+ new = add_refl(out, he, ke, le);
+ copy_data(new, refl);
}
}
- return new;
+ return out;
}
int main(int argc, char *argv[])
{
int c;
- double *input_ref;
- double *phases;
- double *esds;
- char *template = NULL;
int config_noise = 0;
int config_poisson = 0;
- int config_nophase = 0;
int config_multi = 0;
char *holo = NULL;
char *mero = NULL;
char *expand = NULL;
+ char *input_file = NULL;
+ char *template = NULL;
char *output = NULL;
- char *input = NULL;
- char *filename = NULL;
- ReflItemList *input_items;
- ReflItemList *write_items;
- UnitCell *cell = NULL;
char *beamfile = NULL;
struct beam_params *beam = NULL;
+ RefList *input;
+ UnitCell *cell = NULL;
/* Long options */
const struct option longopts[] = {
@@ -290,6 +301,7 @@ int main(int argc, char *argv[])
{"intensities", 1, NULL, 'i'},
{"multiplicity", 0, &config_multi, 1},
{"beam", 1, NULL, 'b'},
+ {"pdb", 1, NULL, 'p'},
{0, 0, NULL, 0}
};
@@ -311,7 +323,7 @@ int main(int argc, char *argv[])
break;
case 'i' :
- input = strdup(optarg);
+ input_file = strdup(optarg);
break;
case 'y' :
@@ -330,6 +342,14 @@ int main(int argc, char *argv[])
beamfile = strdup(optarg);
break;
+ case 'p' :
+ cell = load_cell_from_pdb(optarg);
+ if ( cell == NULL ) {
+ ERROR("Failed to get cell from '%s'\n", optarg);
+ return 1;
+ }
+ break;
+
case 0 :
break;
@@ -339,10 +359,6 @@ int main(int argc, char *argv[])
}
- if ( filename == NULL ) {
- filename = strdup("molecule.pdb");
- }
-
if ( (holo != NULL) && (expand != NULL) ) {
ERROR("You cannot 'twin' and 'expand' at the same time.\n");
ERROR("Decide which one you want to do first.\n");
@@ -358,19 +374,14 @@ int main(int argc, char *argv[])
}
}
- cell = load_cell_from_pdb(filename);
- if ( !config_nophase ) {
- phases = new_list_phase();
- } else {
- phases = NULL;
+ if ( cell == NULL ) {
+ ERROR("You need to give a PDB file with the unit cell.\n");
+ return 1;
}
- esds = new_list_sigma();
- input_ref = new_list_intensity();
- input_items = read_reflections(input, input_ref, phases,
- NULL, esds);
+ input = read_reflections(input_file);
free(input);
- if ( check_symmetry(input_items, mero) ) {
+ if ( check_list_symmetry(input, mero) ) {
ERROR("The input reflection list does not appear to"
" have symmetry %s\n", mero);
return 1;
@@ -378,8 +389,7 @@ int main(int argc, char *argv[])
if ( config_poisson ) {
if ( beam != NULL ) {
- poisson_reflections(input_ref, input_items,
- beam->adu_per_photon);
+ poisson_reflections(input, beam->adu_per_photon);
} else {
ERROR("You must give a beam parameters file in order"
" to calculate Poisson noise.\n");
@@ -387,72 +397,69 @@ int main(int argc, char *argv[])
}
}
- if ( config_noise ) noise_reflections(input_ref, input_items);
+ if ( config_noise ) noise_reflections(input);
if ( holo != NULL ) {
- ReflItemList *new;
+ RefList *new;
STATUS("Twinning from %s into %s\n", mero, holo);
- new = twin_reflections(input_ref, input_items,
- holo, mero, esds);
- delete_items(input_items);
- input_items = new;
+ new = twin_reflections(input, holo, mero);
+
+ /* Replace old with new */
+ reflist_free(input);
+ input = new;
+
+ /* The symmetry of the list has changed */
+ free(mero);
+ mero = holo;
}
if ( expand != NULL ) {
- ReflItemList *new;
+ RefList *new;
STATUS("Expanding from %s into %s\n", mero, expand);
- new = expand_reflections(input_ref, input_items, expand, mero);
- delete_items(input_items);
- input_items = new;
+ new = expand_reflections(input, expand, mero);
+
+ /* Replace old with new */
+ reflist_free(input);
+ input = new;
}
if ( config_multi ) {
- int i;
+ Reflection *refl;
+ RefListIterator *iter;
- for ( i=0; i<num_items(input_items); i++ ) {
+ for ( refl = first_refl(input, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
- struct refl_item *it;
double inty;
+ signed int h, k, l;
+
+ get_indices(refl, &h, &k, &l);
+ inty = get_intensity(refl);
- it = get_item(input_items, i);
- inty = lookup_intensity(input_ref, it->h, it->k, it->l);
- inty *= num_equivs(it->h, it->k, it->l, mero);
- set_intensity(input_ref, it->h, it->k, it->l, inty);
- STATUS("%i %i %i %i\n", it->h, it->k, it->l,
- num_equivs(it->h, it->k, it->l, mero));
+ inty *= (double)num_equivs(h, k, l, mero);
+ set_int(refl, inty);
}
}
if ( template ) {
- /* Write out only reflections which are in the template
- * (and which we have in the input) */
- ReflItemList *template_items;
- template_items = read_reflections(template,
- NULL, NULL, NULL, NULL);
- write_items = intersection_items(input_items, template_items);
- delete_items(template_items);
-
- } else {
-
- /* Write out all reflections */
- write_items = new_items();
- /* (quick way of copying a list) */
- union_items(write_items, input_items);
+ RefList *t = read_reflections(template);
+ RefList *new = template_reflections(input, t);
+ reflist_free(input);
+ input = new;
}
- write_reflections(output, write_items, input_ref, esds, phases,
- NULL, cell);
+ write_reflist(output, input, cell);
- delete_items(input_items);
- delete_items(write_items);
+ reflist_free(input);
return 0;
}
diff --git a/src/indexamajig.c b/src/indexamajig.c
index e1e0e8fd..c0757edd 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -31,11 +31,11 @@
#include "peaks.h"
#include "detector.h"
#include "filters.h"
-#include "reflections.h"
#include "thread-pool.h"
#include "beam-parameters.h"
#include "geometry.h"
#include "stream.h"
+#include "reflist-utils.h"
enum {
diff --git a/src/partialator.c b/src/partialator.c
index 88926582..553a216a 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -26,7 +26,6 @@
#include "utils.h"
#include "hdf5-file.h"
#include "symmetry.h"
-#include "reflections.h"
#include "stream.h"
#include "geometry.h"
#include "peaks.h"
@@ -35,6 +34,7 @@
#include "post-refinement.h"
#include "hrs-scaling.h"
#include "reflist.h"
+#include "reflist-utils.h"
static void show_help(const char *s)
@@ -418,7 +418,7 @@ int main(int argc, char *argv[])
}
/* Output results */
- write_reflections(outfile, obs, I_full, NULL, NULL, cts, NULL);
+ //write_reflist(outfile, obs, I_full, NULL, NULL, cts, NULL);
/* Clean up */
for ( i=0; i<n_total_patterns; i++ ) {
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 6d3c07ab..ba760e32 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -29,9 +29,10 @@
#include "hdf5-file.h"
#include "detector.h"
#include "peaks.h"
-#include "reflections.h"
#include "beam-parameters.h"
#include "symmetry.h"
+#include "reflist.h"
+#include "reflist-utils.h"
static void show_help(const char *s)
@@ -412,34 +413,29 @@ int main(int argc, char *argv[])
} else {
- int i;
- ReflItemList *items;
+ RefList *reflections;
+
+ reflections = read_reflections(intfile);
+ free(intfile);
if ( grad == GRADIENT_PHASED ) {
- phases = new_list_phase();
+ phases = phases_from_list(reflections);
} else {
phases = NULL;
}
- intensities = new_list_intensity();
- phases = new_list_phase();
- flags = new_list_flag();
- items = read_reflections(intfile, intensities, phases,
- NULL, NULL);
- free(intfile);
-
- for ( i=0; i<num_items(items); i++ ) {
- struct refl_item *it = get_item(items, i);
- set_flag(flags, it->h, it->k, it->l, 1);
- }
+ intensities = intensities_from_list(reflections);
+ phases = phases_from_list(reflections);
+ flags = flags_from_list(reflections);
/* Check that the intensities have the correct symmetry */
- if ( check_symmetry(items, sym) ) {
+ if ( check_list_symmetry(reflections, sym) ) {
ERROR("The input reflection list does not appear to"
" have symmetry %s\n", sym);
return 1;
}
- delete_items(items);
+ reflist_free(reflections);
+
}
image.det = get_detector_geometry(geometry);
diff --git a/src/povray.c b/src/povray.c
index 808d3e22..8d0ec048 100644
--- a/src/povray.c
+++ b/src/povray.c
@@ -23,13 +23,13 @@
#include "utils.h"
#include "symmetry.h"
#include "render_hkl.h"
+#include "povray.h"
#define MAX_PROC (256)
-int povray_render_animation(UnitCell *cell, double *ref, unsigned int *counts,
- ReflItemList *items, unsigned int nproc,
+int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
const char *sym, int wght, double boost)
{
FILE *fh;
@@ -39,6 +39,8 @@ int povray_render_animation(UnitCell *cell, double *ref, unsigned int *counts,
pid_t pids[MAX_PROC];
float max;
int i;
+ Reflection *refl;
+ RefListIterator *iter;
if ( (nproc > MAX_PROC) || (nproc < 1) ) {
ERROR("Number of processes must be a number between 1 and %i\n",
@@ -165,27 +167,29 @@ int povray_render_animation(UnitCell *cell, double *ref, unsigned int *counts,
fprintf(fh, "}\n");
max = 0.0;
- for ( i=0; i<num_items(items); i++ ) {
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
- struct refl_item *it;
float val;
+ signed int h, k, l;
- it = get_item(items, i);
+ get_indices(refl, &h, &k, &l);
switch ( wght ) {
case WGHT_I :
- val = lookup_intensity(ref, it->h, it->k, it->l);
+ val = get_intensity(refl);
break;
case WGHT_SQRTI :
- val = lookup_intensity(ref, it->h, it->k, it->l);
+ val = get_intensity(refl);
val = (val>0.0) ? sqrt(val) : 0.0;
break;
case WGHT_COUNTS :
- val = lookup_count(counts, it->h, it->k, it->l);
- val /= (float)num_equivs(it->h, it->k, it->l, sym);
+ val = get_redundancy(refl);
+ val /= (float)num_equivs(h, k, l, sym);
break;
case WGHT_RAWCOUNTS :
- val = lookup_count(counts, it->h, it->k, it->l);
+ val = get_redundancy(refl);
break;
default :
ERROR("Invalid weighting.\n");
@@ -197,30 +201,31 @@ int povray_render_animation(UnitCell *cell, double *ref, unsigned int *counts,
}
max /= boost;
- for ( i=0; i<num_items(items); i++ ) {
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
- struct refl_item *it;
- float radius;
+ signed int h, k, l;float radius;
int s;
float val, p, r, g, b, trans;
int j;
- it = get_item(items, i);
+ get_indices(refl, &h, &k, &l);
switch ( wght ) {
case WGHT_I :
- val = lookup_intensity(ref, it->h, it->k, it->l);
+ val = get_intensity(refl);
break;
case WGHT_SQRTI :
- val = lookup_intensity(ref, it->h, it->k, it->l);
+ val = get_intensity(refl);
val = (val>0.0) ? sqrt(val) : 0.0;
break;
case WGHT_COUNTS :
- val = lookup_count(counts, it->h, it->k, it->l);
- val /= (float)num_equivs(it->h, it->k, it->l, sym);
+ val = get_redundancy(refl);
+ val /= (float)num_equivs(h, k, l, sym);
break;
case WGHT_RAWCOUNTS :
- val = lookup_count(counts, it->h, it->k, it->l);
+ val = get_redundancy(refl);
break;
default :
ERROR("Invalid weighting.\n");
@@ -269,12 +274,12 @@ int povray_render_animation(UnitCell *cell, double *ref, unsigned int *counts,
trans = 1.0-(val/max);
/* For each equivalent */
- for ( j=0; j<num_equivs(it->h, it->k, it->l, sym); j++ ) {
+ for ( j=0; j<num_equivs(h, k, l, sym); j++ ) {
signed int he, ke, le;
float x, y, z;
- get_equiv(it->h, it->k, it->l, &he, &ke, &le, sym, j);
+ get_equiv(h, k, l, &he, &ke, &le, sym, j);
x = asx*he + bsx*ke + csx*le;
y = asy*he + bsy*ke + csy*le;
diff --git a/src/povray.h b/src/povray.h
index 69bf657e..0b0b3446 100644
--- a/src/povray.h
+++ b/src/povray.h
@@ -3,7 +3,7 @@
*
* Invoke POV-ray
*
- * (c) 2006-2010 Thomas White <taw@physics.org>
+ * (c) 2006-2011 Thomas White <taw@physics.org>
*
* Part of CrystFEL - crystallography with a FEL
*
@@ -16,8 +16,10 @@
#ifndef POVRAY_H
#define POVRAY_H
-extern int povray_render_animation(UnitCell *cell, double *ref,
- unsigned int *counts, ReflItemList *items,
+#include "reflist.h"
+#include "cell.h"
+
+extern int povray_render_animation(UnitCell *cell, RefList *list,
unsigned int nproc, const char *sym,
int wght, double boost);
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 6ff40159..f5e705bd 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -468,11 +468,11 @@ int main(int argc, char *argv[])
output = strdup("processed.hkl");
}
- if ( pdb == NULL ) {
- pdb = strdup("molecule.pdb");
- }
-
cell = load_cell_from_pdb(pdb);
+ if ( cell == NULL ) {
+ ERROR("Failed to load cell from '%s'\n", pdb);
+ return 1;
+ }
free(pdb);
if ( sym == NULL ) sym = strdup("1");
@@ -537,7 +537,6 @@ int main(int argc, char *argv[])
write_reflist(output, model, cell);
- cell_free(cell);
fclose(fh);
free(sym);
reflist_free(model);
diff --git a/src/reflections.c b/src/reflections.c
deleted file mode 100644
index 4743aeee..00000000
--- a/src/reflections.c
+++ /dev/null
@@ -1,230 +0,0 @@
-/*
- * reflections.c
- *
- * Utilities for handling reflections
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-
-#include <stdlib.h>
-#include <math.h>
-#include <stdio.h>
-#include <complex.h>
-#include <string.h>
-
-#include "utils.h"
-#include "cell.h"
-#include "reflections.h"
-#include "beam-parameters.h"
-
-
-double *poisson_esds(double *intensities, ReflItemList *items,
- double adu_per_photon)
-{
- int i;
- double *esds = new_list_sigma();
-
- for ( i=0; i<num_items(items); i++ ) {
-
- struct refl_item *it;
- signed int h, k, l;
- double intensity, sigma;
-
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
-
- intensity = lookup_intensity(intensities, h, k, l);
-
- if ( intensity > 0.0 ) {
- sigma = adu_per_photon * sqrt(intensity/adu_per_photon);
- } else {
- sigma = 0.0;
- }
-
- set_sigma(esds, h, k, l, sigma);
-
- }
-
- return esds;
-}
-
-
-void write_reflections(const char *filename, ReflItemList *items,
- double *intensities, double *esds, double *phases,
- unsigned int *counts, UnitCell *cell)
-{
- FILE *fh;
- int i;
-
- if ( filename == NULL ) {
- fh = stdout;
- } else {
- fh = fopen(filename, "w");
- }
-
- if ( fh == NULL ) {
- ERROR("Couldn't open output file '%s'.\n", filename);
- return;
- }
-
- /* Write spacings and angle if zone axis pattern */
- fprintf(fh, " h k l I phase sigma(I) "
- " 1/d(nm^-1) counts\n");
-
- for ( i=0; i<num_items(items); i++ ) {
-
- struct refl_item *it;
- signed int h, k, l;
-
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
-
- int N;
- double intensity, s, sigma;
- char ph[32];
-
- if ( counts != NULL ) {
- N = lookup_count(counts, h, k, l);
- } else {
- N = 1;
- }
-
- intensity = lookup_intensity(intensities, h, k, l);
-
- if ( phases != NULL ) {
- double p;
- p = lookup_phase(phases, h, k, l);
- snprintf(ph, 31, "%8.6f", p);
- } else {
- strncpy(ph, " -", 31);
- }
-
- if ( cell != NULL ) {
- s = 2.0*resolution(cell, h, k, l);
- } else {
- s = 0.0;
- }
-
- if ( esds != NULL ) {
- sigma = lookup_sigma(esds, h, k, l);
- } else {
- sigma = 0.0;
- }
-
- /* h, k, l, I, sigma(I), s */
- fprintf(fh, "%3i %3i %3i %10.2f %s %10.2f %10.2f %7i\n",
- h, k, l, intensity, ph, sigma, s/1.0e9, N);
-
- }
- fclose(fh);
-}
-
-
-/* Read reflections from file. Returns the list of reflections successfully
- * read in. "intensities", "phases" and "counts" are lists which will be
- * populated with the values read from the file. Existing values in either list
- * will be overwritten if the reflection is read from the file, but other values
- * will be left intact.
- *
- * "intensities", "phases" or "counts" can be NULL, if you don't need them.
- */
-ReflItemList *read_reflections(const char *filename,
- double *intensities, double *phases,
- unsigned int *counts, double *esds)
-{
- FILE *fh;
- char *rval;
- ReflItemList *items;
-
- fh = fopen(filename, "r");
- if ( fh == NULL ) {
- ERROR("Failed to open input file '%s'\n", filename);
- return NULL;
- }
-
- items = new_items();
-
- do {
-
- char line[1024];
- signed int h, k, l;
- float intensity, ph, res, sigma;
- char phs[1024];
- int r;
- int cts;
-
- rval = fgets(line, 1023, fh);
- if ( rval == NULL ) continue;
- r = sscanf(line, "%i %i %i %f %s %f %f %i",
- &h, &k, &l, &intensity, phs, &sigma, &res, &cts);
- if ( r >= 8 ) {
- /* Woohoo */
- } else if ( r >= 7 ) {
- /* No "counts", that's fine.. */
- cts = 1;
- } else if ( r >= 6 ) {
- /* No resolution. Didn't want it anyway. */
- res = 0.0;
- } else if ( r >= 5 ) {
- /* No sigma. It's OK today, but one
- * day I'll get you... */
- sigma = 0.0;
- } else if ( r >= 4 ) {
- /* No phase. Better not need it.. */
- if ( phases != NULL ) {
- ERROR("Need phases and none were specified!\n");
- abort();
- }
- } else {
- /* You lose. */
- continue;
- }
-
- add_item(items, h, k, l);
-
- if ( intensities != NULL ) {
- set_intensity(intensities, h, k, l, intensity);
- }
- if ( phases != NULL ) {
- ph = atof(phs);
- set_phase(phases, h, k, l, ph);
- }
- if ( counts != NULL ) {
- set_count(counts, h, k, l, cts);
- }
- if ( esds != NULL ) {
- set_sigma(esds, h, k, l, sigma);
- }
-
- } while ( rval != NULL );
-
- fclose(fh);
-
- return items;
-}
-
-
-double *ideal_intensities(double complex *sfac)
-{
- double *ref;
- signed int h, k, l;
-
- ref = new_list_intensity();
-
- /* Generate ideal reflections from complex structure factors */
- for ( h=-INDMAX; h<=INDMAX; h++ ) {
- for ( k=-INDMAX; k<=INDMAX; k++ ) {
- for ( l=-INDMAX; l<=INDMAX; l++ ) {
- double complex F = lookup_sfac(sfac, h, k, l);
- double intensity = pow(cabs(F), 2.0);
- set_intensity(ref, h, k, l, intensity);
- }
- }
- }
-
- return ref;
-}
diff --git a/src/reflections.h b/src/reflections.h
deleted file mode 100644
index 20060f21..00000000
--- a/src/reflections.h
+++ /dev/null
@@ -1,38 +0,0 @@
-/*
- * reflections.h
- *
- * Utilities for handling reflections
- *
- * (c) 2006-2010 Thomas White <taw@physics.org>
- *
- * Part of CrystFEL - crystallography with a FEL
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-#ifndef REFLECTIONS_H
-#define REFLECTIONS_H
-
-
-#include "cell.h"
-#include "utils.h"
-
-
-extern double *poisson_esds(double *intensities, ReflItemList *items,
- double adu_per_photon);
-
-extern void write_reflections(const char *filename, ReflItemList *items,
- double *intensities, double *esds, double *phases,
- unsigned int *counts, UnitCell *cell);
-
-extern ReflItemList *read_reflections(const char *filename,
- double *intensities, double *phases,
- unsigned int *counts, double *esds);
-
-extern double *ideal_intensities(double complex *sfac);
-
-
-#endif /* REFLECTIONS_H */
diff --git a/src/reflist-utils.c b/src/reflist-utils.c
index 33d373b8..aac94fbf 100644
--- a/src/reflist-utils.c
+++ b/src/reflist-utils.c
@@ -15,6 +15,153 @@
#include "reflist.h"
#include "cell.h"
+#include "utils.h"
+#include "reflist-utils.h"
+#include "symmetry.h"
+
+
+double *intensities_from_list(RefList *list)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+ double *out = new_list_intensity();
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
+
+ signed int h, k, l;
+ double intensity = get_intensity(refl);
+
+ get_indices(refl, &h, &k, &l);
+
+ set_intensity(out, h, k, l, intensity);
+
+ }
+
+ return out;
+}
+
+
+double *phases_from_list(RefList *list)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+ double *out = new_list_phase();
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
+
+ signed int h, k, l;
+ double phase = get_phase(refl);
+
+ get_indices(refl, &h, &k, &l);
+
+ set_phase(out, h, k, l, phase);
+
+ }
+
+ return out;
+
+}
+
+
+unsigned char *flags_from_list(RefList *list)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+ unsigned char *out = new_list_flag();
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
+
+ signed int h, k, l;
+
+ get_indices(refl, &h, &k, &l);
+
+ set_flag(out, h, k, l, 1);
+
+ }
+
+ return out;
+
+}
+
+
+int check_list_symmetry(RefList *list, const char *sym)
+{
+ unsigned char *flags;
+ Reflection *refl;
+ RefListIterator *iter;
+
+ flags = flags_from_list(list);
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
+
+ int j;
+ int found = 0;
+ signed int h, k, l;
+
+ get_indices(refl, &h, &k, &l);
+
+ for ( j=0; j<num_equivs(h, k, l, sym); j++ ) {
+
+ signed int he, ke, le;
+ get_equiv(h, k, l, &he, &ke, &le, sym, j);
+
+ if ( abs(he) > INDMAX ) continue;
+ if ( abs(le) > INDMAX ) continue;
+ if ( abs(ke) > INDMAX ) continue;
+
+ found += lookup_flag(flags, he, ke, le);
+
+ }
+
+ if ( found > 1 ) {
+ free(flags);
+ return 1; /* Symmetry is wrong! */
+ }
+
+ }
+
+ free(flags);
+
+ return 0;
+}
+
+
+int find_equiv_in_list(RefList *list, signed int h, signed int k,
+ signed int l, const char *sym, signed int *hu,
+ signed int *ku, signed int *lu)
+{
+ int i;
+ int found = 0;
+
+ for ( i=0; i<num_equivs(h, k, l, sym); i++ ) {
+
+ signed int he, ke, le;
+ Reflection *f;
+ get_equiv(h, k, l, &he, &ke, &le, sym, i);
+ f = find_refl(list, he, ke, le);
+
+ /* There must only be one equivalent. If there are more, it
+ * indicates that the user lied about the input symmetry.
+ * This situation should have been checked for earlier by
+ * calling check_symmetry() with 'items' and 'mero'. */
+
+ if ( (f != NULL) && !found ) {
+ *hu = he; *ku = ke; *lu = le;
+ return 1;
+ }
+
+ }
+
+ return 0;
+}
void write_reflections_to_file(FILE *fh, RefList *list, UnitCell *cell)
@@ -31,18 +178,26 @@ void write_reflections_to_file(FILE *fh, RefList *list, UnitCell *cell)
signed int h, k, l;
double intensity, esd_i, s;
+ int red;
double fs, ss;
+ char res[16];
get_indices(refl, &h, &k, &l);
get_detector_pos(refl, &fs, &ss);
intensity = get_intensity(refl);
- esd_i = 0.0; /* FIXME! */
- s = 0.0; /* FIXME! */
+ esd_i = get_esd_intensity(refl);
+ red = get_redundancy(refl);
+
+ if ( cell != NULL ) {
+ s = resolution(cell, h, k, l);
+ snprintf(res, 16, "%10.2f", s/1e9);
+ } else {
+ strcpy(res, " -");
+ }
- /* h, k, l, I, sigma(I), s */
fprintf(fh,
- "%3i %3i %3i %10.2f %s %10.2f %10.2f %7i %6.1f %6.1f\n",
- h, k, l, intensity, " -", esd_i, s/1.0e9, 1,
+ "%3i %3i %3i %10.2f %s %10.2f %s %7i %6.1f %6.1f\n",
+ h, k, l, intensity, " -", esd_i, res, red,
fs, ss);
}
@@ -70,3 +225,85 @@ int write_reflist(const char *filename, RefList *list, UnitCell *cell)
return 0;
}
+
+
+RefList *read_reflections_from_file(FILE *fh)
+{
+ char *rval = NULL;
+ int first = 1;
+ RefList *out;
+
+ out = reflist_new();
+
+ do {
+
+ char line[1024];
+ signed int h, k, l;
+ float intensity, sigma, fs, ss;
+ char phs[1024];
+ char ress[1024];
+ int cts;
+ int r;
+ Reflection *refl;
+
+ rval = fgets(line, 1023, fh);
+ if ( rval == NULL ) continue;
+ chomp(line);
+
+ if ( strcmp(line, REFLECTION_END_MARKER) == 0 ) return out;
+
+ r = sscanf(line, "%i %i %i %f %s %f %s %i %f %f",
+ &h, &k, &l, &intensity, phs, &sigma, ress, &cts,
+ &fs, &ss);
+ if ( (r != 10) && (!first) ) {
+ reflist_free(out);
+ return NULL;
+ }
+
+ first = 0;
+ if ( r == 10 ) {
+
+ double ph;
+ char *v;
+
+ refl = add_refl(out, h, k, l);
+ set_int(refl, intensity);
+ set_detector_pos(refl, fs, ss, 0.0);
+ set_esd_intensity(refl, sigma);
+
+ ph = strtod(phs, &v);
+ if ( v != NULL ) set_ph(refl, ph);
+
+ /* The 1/d value is actually ignored. */
+
+ }
+
+ } while ( rval != NULL );
+
+ /* Got read error of some kind before finding PEAK_LIST_END_MARKER */
+ return NULL;
+}
+
+
+RefList *read_reflections(const char *filename)
+{
+ FILE *fh;
+ RefList *out;
+
+ if ( filename == NULL ) {
+ fh = stdout;
+ } else {
+ fh = fopen(filename, "r");
+ }
+
+ if ( fh == NULL ) {
+ ERROR("Couldn't open input file '%s'.\n", filename);
+ return NULL;
+ }
+
+ out = read_reflections_from_file(fh);
+
+ fclose(fh);
+
+ return out;
+}
diff --git a/src/reflist-utils.h b/src/reflist-utils.h
index d5c1e655..b7ecd3ea 100644
--- a/src/reflist-utils.h
+++ b/src/reflist-utils.h
@@ -21,10 +21,24 @@
#include "cell.h"
+#define REFLECTION_END_MARKER "End of reflections"
+
+
extern void write_reflections_to_file(FILE *fh, RefList *list, UnitCell *cell);
-extern int write_reflist(const char *filename, RefList *list,
- UnitCell *cell);
+extern int write_reflist(const char *filename, RefList *list, UnitCell *cell);
+
+extern RefList *read_reflections_from_file(FILE *fh);
+
+extern RefList *read_reflections(const char *filename);
+
+extern double *intensities_from_list(RefList *list);
+extern double *phases_from_list(RefList *list);
+extern unsigned char *flags_from_list(RefList *list);
+extern int check_list_symmetry(RefList *list, const char *sym);
+extern int find_equiv_in_list(RefList *list, signed int h, signed int k,
+ signed int l, const char *sym, signed int *hu,
+ signed int *ku, signed int *lu);
#endif /* REFLIST_UTILS_H */
diff --git a/src/reflist.c b/src/reflist.c
index 3181b418..a81c738e 100644
--- a/src/reflist.c
+++ b/src/reflist.c
@@ -46,6 +46,9 @@ struct _refldata {
double intensity;
double esd_i;
+ /* Phase */
+ double phase;
+
/* Redundancy */
int redundancy;
@@ -250,8 +253,20 @@ double get_esd_intensity(Reflection *refl)
}
+double get_phase(Reflection *refl)
+{
+ return refl->data.phase;
+}
+
+
/********************************** Setters ***********************************/
+void copy_data(Reflection *to, Reflection *from)
+{
+ memcpy(&to->data, &from->data, sizeof(struct _refldata));
+}
+
+
void set_detector_pos(Reflection *refl, double exerr, double x, double y)
{
refl->data.excitation_error = exerr;
@@ -301,6 +316,12 @@ void set_esd_intensity(Reflection *refl, double esd)
}
+void set_ph(Reflection *refl, double phase)
+{
+ refl->data.phase = phase;
+}
+
+
/********************************* Insertion **********************************/
static void insert_node(Reflection *head, Reflection *new)
@@ -689,3 +710,9 @@ void optimise_reflist(RefList *list)
recursive_depth(list->head->child[0]));
}
}
+
+
+int num_reflections(RefList *list)
+{
+ return recursive_count(list->head->child[0]);
+}
diff --git a/src/reflist.h b/src/reflist.h
index f6d01ecf..5b647102 100644
--- a/src/reflist.h
+++ b/src/reflist.h
@@ -44,8 +44,10 @@ extern int get_scalable(Reflection *refl);
extern int get_redundancy(Reflection *refl);
extern double get_sum_squared_dev(Reflection *refl);
extern double get_esd_intensity(Reflection *refl);
+extern double get_phase(Reflection *refl);
/* Set */
+extern void copy_data(Reflection *to, Reflection *from);
extern void set_detector_pos(Reflection *refl, double exerr,
double x, double y);
extern void set_partial(Reflection *refl, double r1, double r2, double p,
@@ -55,6 +57,7 @@ extern void set_scalable(Reflection *refl, int scalable);
extern void set_redundancy(Reflection *refl, int red);
extern void set_sum_squared_dev(Reflection *refl, double dev);
extern void set_esd_intensity(Reflection *refl, double esd);
+extern void set_ph(Reflection *refl, double phase);
/* Insertion */
extern Reflection *add_refl(RefList *list, INDICES);
@@ -66,7 +69,8 @@ extern void delete_refl(Reflection *refl);
extern Reflection *first_refl(RefList *list, RefListIterator **iterator);
extern Reflection *next_refl(Reflection *refl, RefListIterator *iter);
-/* Voodoo */
+/* Misc */
extern void optimise_reflist(RefList *list);
+extern int num_reflections(RefList *list);
#endif /* REFLIST_H */
diff --git a/src/render_hkl.c b/src/render_hkl.c
index fc1c5121..aec6d714 100644
--- a/src/render_hkl.c
+++ b/src/render_hkl.c
@@ -26,11 +26,12 @@
#endif
#include "utils.h"
-#include "reflections.h"
#include "povray.h"
#include "symmetry.h"
#include "render.h"
#include "render_hkl.h"
+#include "reflist.h"
+#include "reflist-utils.h"
static void show_help(const char *s)
@@ -81,8 +82,7 @@ static void show_help(const char *s)
static void draw_circles(signed int xh, signed int xk, signed int xl,
signed int yh, signed int yk, signed int yl,
signed int zh, signed int zk, signed int zl,
- double *ref, unsigned int *counts, ReflItemList *items,
- const char *sym,
+ RefList *list, const char *sym,
cairo_t *dctx, int wght, double boost, int colscale,
UnitCell *cell, double radius, double theta,
double as, double bs, double cx, double cy,
@@ -98,36 +98,36 @@ static void draw_circles(signed int xh, signed int xk, signed int xl,
*max_res = 0.0; *max_ux = 0; *max_uy = 0;
}
- /* Loop across the two basis directions */
+ /* Iterate over possible reflections in this zone */
for ( xi=-INDMAX; xi<INDMAX; xi++ ) {
for ( yi=-INDMAX; yi<INDMAX; yi++ ) {
double u, v, val, res;
signed int h, k, l;
- signed int he, ke, le;
+ Reflection *refl;
h = xi*xh + yi*yh;
k = xi*xk + yi*yk;
l = xi*xl + yi*yl;
/* Got this reflection? */
- if ( find_unique_equiv(items, h, k, l, sym,
- &he, &ke, &le) == 0 ) continue;
+ refl = find_refl(list, h, k, l);
+ if ( refl == NULL ) continue;
switch ( wght ) {
case WGHT_I :
- val = lookup_intensity(ref, he, ke, le);
+ val = get_intensity(refl);
break;
case WGHT_SQRTI :
- val = lookup_intensity(ref, he, ke, le);
+ val = get_intensity(refl);
val = (val>0.0) ? sqrt(val) : 0.0;
break;
case WGHT_COUNTS :
- val = lookup_count(counts, he, ke, le);
- val /= (float)num_equivs(he, ke, le, sym);
+ val = get_redundancy(refl);
+ val /= (float)num_equivs(h, k, l, sym);
break;
case WGHT_RAWCOUNTS :
- val = lookup_count(counts, he, ke, le);
+ val = get_redundancy(refl);
break;
default :
ERROR("Invalid weighting.\n");
@@ -228,8 +228,7 @@ static void render_overlined_indices(cairo_t *dctx,
}
-static void render_za(UnitCell *cell, ReflItemList *items,
- double *ref, unsigned int *counts,
+static void render_za(UnitCell *cell, RefList *list,
double boost, const char *sym, int wght, int colscale,
signed int xh, signed int xk, signed int xl,
signed int yh, signed int yk, signed int yl,
@@ -284,7 +283,7 @@ static void render_za(UnitCell *cell, ReflItemList *items,
scale = 1.0;
draw_circles(xh, xk, xl, yh, yk, yl, zh, zk, zl,
- ref, counts, items, sym, NULL, wght, boost, colscale, cell,
+ list, sym, NULL, wght, boost, colscale, cell,
0.0, theta, as, bs, 0.0, 0.0, scale,
&max_ux, &max_uy, &max_val, &max_u, &max_v, &max_res);
@@ -342,7 +341,7 @@ static void render_za(UnitCell *cell, ReflItemList *items,
cy = 512.0 - 20.0;
draw_circles(xh, xk, xl, yh, yk, yl, zh, zk, zl,
- ref, counts, items, sym, dctx, wght, boost, colscale, cell,
+ list, sym, dctx, wght, boost, colscale, cell,
max_r, theta, as, bs, cx, cy, scale,
NULL, NULL, &max_val, NULL, NULL, NULL);
@@ -462,8 +461,8 @@ int main(int argc, char *argv[])
{
int c;
UnitCell *cell;
+ RefList *list;
char *infile;
- double *ref;
int config_povray = 0;
int config_zoneaxis = 0;
int config_sqrt = 0;
@@ -477,7 +476,6 @@ int main(int argc, char *argv[])
int wght;
int colscale;
char *cscale = NULL;
- unsigned int *cts;
signed int dh=1, dk=0, dl=0;
signed int rh=0, rk=1, rl=0;
char *down = NULL;
@@ -613,7 +611,8 @@ int main(int argc, char *argv[])
if ( config_zoneaxis ) {
if ( (( down == NULL ) && ( right != NULL ))
|| (( down != NULL ) && ( right == NULL )) ) {
- ERROR("Either specify both 'down' and 'right', or neither.\n");
+ ERROR("Either specify both 'down' and 'right',"
+ " or neither.\n");
return 1;
}
if ( down != NULL ) {
@@ -641,24 +640,22 @@ int main(int argc, char *argv[])
ERROR("Couldn't load unit cell from %s\n", pdb);
return 1;
}
- ref = new_list_intensity();
- cts = new_list_count();
- ReflItemList *items = read_reflections(infile, ref, NULL, cts, NULL);
- if ( items == NULL ) {
- ERROR("Couldn't open file '%s'\n", infile);
+ list = read_reflections(infile);
+ if ( list == NULL ) {
+ ERROR("Couldn't read file '%s'\n", infile);
return 1;
}
- if ( check_symmetry(items, sym) ) {
+ if ( check_list_symmetry(list, sym) ) {
ERROR("The input reflection list does not appear to"
" have symmetry %s\n", sym);
return 1;
}
if ( config_povray ) {
- r = povray_render_animation(cell, ref, cts, items,
+ r = povray_render_animation(cell, list,
nproc, sym, wght, boost);
} else if ( config_zoneaxis ) {
- render_za(cell, items, ref, cts, boost, sym, wght, colscale,
+ render_za(cell, list, boost, sym, wght, colscale,
rh, rk, rl, dh, dk, dl, outfile);
} else {
ERROR("Try again with either --povray or --zone-axis.\n");
@@ -666,7 +663,7 @@ int main(int argc, char *argv[])
free(pdb);
free(sym);
- delete_items(items);
+ reflist_free(list);
if ( outfile != NULL ) free(outfile);
return r;
diff --git a/src/statistics.c b/src/statistics.c
index d0d3ae85..77e7919f 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -25,9 +25,8 @@
struct r_params {
- const double *ref1;
- const double *ref2;
- ReflItemList *items; /* Which reflections to use */
+ RefList *list1;
+ RefList *list2;
int fom; /* Which FoM to use (see the enum just below) */
};
@@ -43,27 +42,30 @@ enum {
/* Return the least squares optimal scaling factor when comparing intensities.
- * ref1,ref2 are the two intensity lists to compare. "items" is a ReflItemList
+ * list1,list2 are the two intensity lists to compare. "items" is a ReflItemList
* containing the reflections which should be taken into account.
*/
-double stat_scale_intensity(const double *ref1, const double *ref2,
- ReflItemList *items)
+double stat_scale_intensity(RefList *list1, RefList *list2)
{
double top = 0.0;
double bot = 0.0;
- int i;
+ Reflection *refl1;
+ RefListIterator *iter;
- for ( i=0; i<num_items(items); i++ ) {
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) ) {
double i1, i2;
- struct refl_item *it;
signed int h, k, l;
+ Reflection *refl2;
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
- i1 = lookup_intensity(ref1, h, k, l);
- i2 = lookup_intensity(ref2, h, k, l);
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
top += i1 * i2;
bot += i2 * i2;
@@ -77,30 +79,36 @@ double stat_scale_intensity(const double *ref1, const double *ref2,
/* Return the least squares optimal scaling factor when comparing the square
* roots of the intensities (i.e. one approximation to the structure factor
* moduli).
- * ref1,ref2 are the two intensity lists to compare (they contain intensities,
+ * list1,list2 are the two intensity lists to compare (they contain intensities,
* not square rooted intensities). "items" is a ReflItemList containing the
* reflections which should be taken into account.
*/
-static double stat_scale_sqrti(const double *ref1, const double *ref2,
- ReflItemList *items)
+static double stat_scale_sqrti(RefList *list1, RefList *list2)
{
double top = 0.0;
double bot = 0.0;
- int i;
+ Reflection *refl1;
+ RefListIterator *iter;
- for ( i=0; i<num_items(items); i++ ) {
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) ) {
- double i1, i2, f1, f2;
- struct refl_item *it;
+ double i1, i2;
+ double f1, f2;
signed int h, k, l;
+ Reflection *refl2;
+
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
- i1 = lookup_intensity(ref1, h, k, l);
if ( i1 < 0.0 ) continue;
f1 = sqrt(i1);
- i2 = lookup_intensity(ref2, h, k, l);
+
if ( i2 < 0.0 ) continue;
f2 = sqrt(i2);
@@ -113,26 +121,33 @@ static double stat_scale_sqrti(const double *ref1, const double *ref2,
}
-static double internal_r1_ignorenegs(const double *ref1, const double *ref2,
- ReflItemList *items, double scale)
+static double internal_r1_ignorenegs(RefList *list1, RefList *list2,
+ double scale)
{
double top = 0.0;
double bot = 0.0;
- int i;
+ Reflection *refl1;
+ RefListIterator *iter;
- for ( i=0; i<num_items(items); i++ ) {
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) ) {
- double i1, f1, i2, f2;
- struct refl_item *it;
+ double i1, i2;
+ double f1, f2;
signed int h, k, l;
+ Reflection *refl2;
+
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
- i1 = lookup_intensity(ref1, h, k, l);
if ( i1 < 0.0 ) continue;
f1 = sqrt(i1);
- i2 = lookup_intensity(ref2, h, k, l);
+
if ( i2 < 0.0 ) continue;
f2 = sqrt(i2);
f2 *= scale;
@@ -146,25 +161,32 @@ static double internal_r1_ignorenegs(const double *ref1, const double *ref2,
}
-static double internal_r1_negstozero(const double *ref1, const double *ref2,
- ReflItemList *items, double scale)
+static double internal_r1_negstozero(RefList *list1, RefList *list2,
+ double scale)
{
double top = 0.0;
double bot = 0.0;
- int i;
+ Reflection *refl1;
+ RefListIterator *iter;
- for ( i=0; i<num_items(items); i++ ) {
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) ) {
- double i1, f1, i2, f2;
- struct refl_item *it;
+ double i1, i2;
+ double f1, f2;
signed int h, k, l;
+ Reflection *refl2;
+
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
- i1 = lookup_intensity(ref1, h, k, l);
f1 = i1 > 0.0 ? sqrt(i1) : 0.0;
- i2 = lookup_intensity(ref2, h, k, l);
+
f2 = i2 > 0.0 ? sqrt(i2) : 0.0;
f2 *= scale;
@@ -177,24 +199,27 @@ static double internal_r1_negstozero(const double *ref1, const double *ref2,
}
-static double internal_r2(const double *ref1, const double *ref2,
- ReflItemList *items, double scale)
+static double internal_r2(RefList *list1, RefList *list2, double scale)
{
double top = 0.0;
double bot = 0.0;
- int i;
+ Reflection *refl1;
+ RefListIterator *iter;
- for ( i=0; i<num_items(items); i++ ) {
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) ) {
double i1, i2;
- struct refl_item *it;
signed int h, k, l;
+ Reflection *refl2;
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
- i1 = lookup_intensity(ref1, h, k, l);
- i2 = scale * lookup_intensity(ref2, h, k, l);
+ i1 = get_intensity(refl1);
+ i2 = scale * get_intensity(refl2);
top += pow(i1 - i2, 2.0);
bot += pow(i1, 2.0);
@@ -205,24 +230,27 @@ static double internal_r2(const double *ref1, const double *ref2,
}
-static double internal_r_i(const double *ref1, const double *ref2,
- ReflItemList *items, double scale)
+static double internal_r_i(RefList *list1, RefList *list2, double scale)
{
double top = 0.0;
double bot = 0.0;
- int i;
+ Reflection *refl1;
+ RefListIterator *iter;
- for ( i=0; i<num_items(items); i++ ) {
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) ) {
double i1, i2;
- struct refl_item *it;
signed int h, k, l;
+ Reflection *refl2;
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
- i1 = lookup_intensity(ref1, h, k, l);
- i2 = scale * lookup_intensity(ref2, h, k, l);
+ i1 = get_intensity(refl1);
+ i2 = scale * get_intensity(refl2);
top += fabs(i1-i2);
bot += fabs(i1);
@@ -233,24 +261,29 @@ static double internal_r_i(const double *ref1, const double *ref2,
}
-static double internal_rdiff_intensity(const double *ref1, const double *ref2,
- ReflItemList *items, double scale)
+static double internal_rdiff_intensity(RefList *list1, RefList *list2,
+ double scale)
{
double top = 0.0;
double bot = 0.0;
- int i;
+ Reflection *refl1;
+ RefListIterator *iter;
- for ( i=0; i<num_items(items); i++ ) {
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) ) {
double i1, i2;
- struct refl_item *it;
signed int h, k, l;
+ Reflection *refl2;
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
+
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
- i1 = lookup_intensity(ref1, h, k, l);
- i2 = lookup_intensity(ref2, h, k, l);
i2 *= scale;
top += fabs(i1 - i2);
@@ -262,25 +295,32 @@ static double internal_rdiff_intensity(const double *ref1, const double *ref2,
}
-static double internal_rdiff_negstozero(const double *ref1, const double *ref2,
- ReflItemList *items, double scale)
+static double internal_rdiff_negstozero(RefList *list1, RefList *list2,
+ double scale)
{
double top = 0.0;
double bot = 0.0;
- int i;
+ Reflection *refl1;
+ RefListIterator *iter;
- for ( i=0; i<num_items(items); i++ ) {
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) ) {
- double i1, i2, f1, f2;
- struct refl_item *it;
+ double i1, i2;
+ double f1, f2;
signed int h, k, l;
+ Reflection *refl2;
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
+
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
- i1 = lookup_intensity(ref1, h, k, l);
f1 = i1 > 0.0 ? sqrt(i1) : 0.0;
- i2 = lookup_intensity(ref2, h, k, l);
+
f2 = i2 > 0.0 ? sqrt(i2) : 0.0;
f2 *= scale;
@@ -293,26 +333,33 @@ static double internal_rdiff_negstozero(const double *ref1, const double *ref2,
}
-static double internal_rdiff_ignorenegs(const double *ref1, const double *ref2,
- ReflItemList *items, double scale)
+static double internal_rdiff_ignorenegs(RefList *list1, RefList *list2,
+ double scale)
{
double top = 0.0;
double bot = 0.0;
- int i;
+ Reflection *refl1;
+ RefListIterator *iter;
- for ( i=0; i<num_items(items); i++ ) {
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) ) {
- double i1, i2, f1, f2;
- struct refl_item *it;
+ double i1, i2;
+ double f1, f2;
signed int h, k, l;
+ Reflection *refl2;
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
+
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
- i1 = lookup_intensity(ref1, h, k, l);
if ( i1 < 0.0 ) continue;
f1 = sqrt(i1);
- i2 = lookup_intensity(ref2, h, k, l);
+
if ( i2 < 0.0 ) continue;
f2 = sqrt(i2);
f2 *= scale;
@@ -332,26 +379,21 @@ static double calc_r(double scale, void *params)
switch ( rp->fom) {
case R_1_ZERO :
- return internal_r1_negstozero(rp->ref1, rp->ref2,
- rp->items, scale);
+ return internal_r1_negstozero(rp->list1, rp->list2, scale);
case R_1_IGNORE :
- return internal_r1_ignorenegs(rp->ref1, rp->ref2,
- rp->items, scale);
+ return internal_r1_ignorenegs(rp->list1, rp->list2, scale);
case R_2 :
- return internal_r2(rp->ref1, rp->ref2, rp->items, scale);
+ return internal_r2(rp->list1, rp->list2, scale);
case R_1_I :
- return internal_r_i(rp->ref1, rp->ref2, rp->items, scale);
+ return internal_r_i(rp->list1, rp->list2, scale);
case R_DIFF_ZERO :
- return internal_rdiff_negstozero(rp->ref1, rp->ref2,
- rp->items, scale);
+ return internal_rdiff_negstozero(rp->list1, rp->list2,scale);
case R_DIFF_IGNORE :
- return internal_rdiff_ignorenegs(rp->ref1, rp->ref2,
- rp->items, scale);
+ return internal_rdiff_ignorenegs(rp->list1, rp->list2, scale);
case R_DIFF_INTENSITY :
- return internal_rdiff_intensity(rp->ref1, rp->ref2,
- rp->items, scale);
+ return internal_rdiff_intensity(rp->list1, rp->list2, scale);
}
ERROR("No such FoM!\n");
@@ -359,8 +401,8 @@ static double calc_r(double scale, void *params)
}
-static double r_minimised(const double *ref1, const double *ref2,
- ReflItemList *items, double *scalep, int fom)
+static double r_minimised(RefList *list1, RefList *list2,
+ double *scalep, int fom)
{
gsl_function F;
gsl_min_fminimizer *s;
@@ -369,9 +411,8 @@ static double r_minimised(const double *ref1, const double *ref2,
struct r_params rp;
int iter = 0;
- rp.ref1 = ref1;
- rp.ref2 = ref2;
- rp.items = items;
+ rp.list1 = list1;
+ rp.list2 = list2;
rp.fom = fom;
F.function = &calc_r;
@@ -385,12 +426,12 @@ static double r_minimised(const double *ref1, const double *ref2,
case R_1_IGNORE :
case R_DIFF_ZERO :
case R_DIFF_IGNORE :
- scale = stat_scale_sqrti(ref1, ref2, items);
+ scale = stat_scale_sqrti(list1, list2);
break;
case R_2 :
case R_1_I :
case R_DIFF_INTENSITY :
- scale = stat_scale_intensity(ref1, ref2, items);
+ scale = stat_scale_intensity(list1, list2);
break;
}
//STATUS("Initial scale factor estimate: %5.2e\n", scale);
@@ -429,77 +470,74 @@ static double r_minimised(const double *ref1, const double *ref2,
}
-double stat_r1_ignore(const double *ref1, const double *ref2,
- ReflItemList *items, double *scalep)
+double stat_r1_ignore(RefList *list1, RefList *list2, double *scalep)
{
- return r_minimised(ref1, ref2, items, scalep, R_1_IGNORE);
+ return r_minimised(list1, list2, scalep, R_1_IGNORE);
}
-double stat_r1_zero(const double *ref1, const double *ref2,
- ReflItemList *items, double *scalep)
+double stat_r1_zero(RefList *list1, RefList *list2, double *scalep)
{
- return r_minimised(ref1, ref2, items, scalep, R_1_ZERO);
+ return r_minimised(list1, list2, scalep, R_1_ZERO);
}
-double stat_r2(const double *ref1, const double *ref2,
- ReflItemList *items, double *scalep)
+double stat_r2(RefList *list1, RefList *list2, double *scalep)
{
- return r_minimised(ref1, ref2, items, scalep, R_2);
+ return r_minimised(list1, list2, scalep, R_2);
}
-double stat_r1_i(const double *ref1, const double *ref2,
- ReflItemList *items, double *scalep)
+double stat_r1_i(RefList *list1, RefList *list2, double *scalep)
{
- return r_minimised(ref1, ref2, items, scalep, R_1_I);
+ return r_minimised(list1, list2, scalep, R_1_I);
}
-double stat_rdiff_zero(const double *ref1, const double *ref2,
- ReflItemList *items, double *scalep)
+double stat_rdiff_zero(RefList *list1, RefList *list2, double *scalep)
{
- return r_minimised(ref1, ref2, items, scalep, R_DIFF_ZERO);
+ return r_minimised(list1, list2, scalep, R_DIFF_ZERO);
}
-double stat_rdiff_ignore(const double *ref1, const double *ref2,
- ReflItemList *items, double *scalep)
+double stat_rdiff_ignore(RefList *list1, RefList *list2, double *scalep)
{
- return r_minimised(ref1, ref2, items, scalep, R_DIFF_IGNORE);
+ return r_minimised(list1, list2, scalep, R_DIFF_IGNORE);
}
-double stat_rdiff_intensity(const double *ref1, const double *ref2,
- ReflItemList *items, double *scalep)
+double stat_rdiff_intensity(RefList *list1, RefList *list2, double *scalep)
{
- return r_minimised(ref1, ref2, items, scalep, R_DIFF_INTENSITY);
+ return r_minimised(list1, list2, scalep, R_DIFF_INTENSITY);
}
-double stat_pearson_i(const double *ref1, const double *ref2,
- ReflItemList *items)
+
+double stat_pearson_i(RefList *list1, RefList *list2)
{
double *vec1, *vec2;
- int i = 0;
- int ni = num_items(items);
+ int ni = num_reflections(list1) + num_reflections(list2);
double val;
int nacc = 0;
+ Reflection *refl1;
+ RefListIterator *iter;
vec1 = malloc(ni*sizeof(double));
vec2 = malloc(ni*sizeof(double));
- for ( i=0; i<ni; i++ ) {
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) ) {
double i1, i2;
- struct refl_item *it;
signed int h, k, l;
+ Reflection *refl2;
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
- i1 = lookup_intensity(ref1, h, k, l);
- i2 = lookup_intensity(ref2, h, k, l);
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
vec1[nacc] = i1;
vec2[nacc] = i2;
@@ -514,40 +552,42 @@ double stat_pearson_i(const double *ref1, const double *ref2,
}
-double stat_pearson_f_ignore(const double *ref1, const double *ref2,
- ReflItemList *items)
+double stat_pearson_f_ignore(RefList *list1, RefList *list2)
{
double *vec1, *vec2;
- int i = 0;
- int ni = num_items(items);
+ int ni = num_reflections(list1) + num_reflections(list2);
double val;
int nacc = 0;
+ Reflection *refl1;
+ RefListIterator *iter;
vec1 = malloc(ni*sizeof(double));
vec2 = malloc(ni*sizeof(double));
- for ( i=0; i<ni; i++ ) {
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) ) {
- double i1, i2, f1, f2;
- struct refl_item *it;
+ double i1, i2;
+ double f1, f2;
signed int h, k, l;
+ Reflection *refl2;
int ig = 0;
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
+
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
- i1 = lookup_intensity(ref1, h, k, l);
- if ( i1 < 0.0 ) {
- ig=1;
- }
+ if ( i1 < 0.0 ) ig = 1;
f1 = sqrt(i1);
- i2 = lookup_intensity(ref2, h, k, l);
- if ( i2 < 0.0 ) {
- ig=1;
- }
+
+ if ( i2 < 0.0 ) ig = 1;
f2 = sqrt(i2);
- if ( ig ) continue;
+ if ( ig ) continue;
vec1[nacc] = f1;
vec2[nacc] = f2;
@@ -563,30 +603,35 @@ double stat_pearson_f_ignore(const double *ref1, const double *ref2,
}
-double stat_pearson_f_zero(const double *ref1, const double *ref2,
- ReflItemList *items)
+double stat_pearson_f_zero(RefList *list1, RefList *list2)
{
double *vec1, *vec2;
- int i = 0;
- int ni = num_items(items);
+ int ni = num_reflections(list1) + num_reflections(list2);
double val;
int nacc = 0;
+ Reflection *refl1;
+ RefListIterator *iter;
vec1 = malloc(ni*sizeof(double));
vec2 = malloc(ni*sizeof(double));
- for ( i=0; i<ni; i++ ) {
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) ) {
- double i1, i2, f1, f2;
- struct refl_item *it;
+ double i1, i2;
+ double f1, f2;
signed int h, k, l;
+ Reflection *refl2;
+
+ get_indices(refl1, &h, &k, &l);
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue; /* No common reflection */
- it = get_item(items, i);
- h = it->h; k = it->k; l = it->l;
+ i1 = get_intensity(refl1);
+ i2 = get_intensity(refl2);
- i1 = lookup_intensity(ref1, h, k, l);
f1 = i1 > 0.0 ? sqrt(i1) : 0.0;
- i2 = lookup_intensity(ref2, h, k, l);
f2 = i2 > 0.0 ? sqrt(i2) : 0.0;
vec1[nacc] = f1;
diff --git a/src/statistics.h b/src/statistics.h
index 18603357..34ef2406 100644
--- a/src/statistics.h
+++ b/src/statistics.h
@@ -17,35 +17,25 @@
#define STATISTICS_H
-#include "utils.h"
-
-extern double stat_scale_intensity(const double *ref1, const double *ref2,
- ReflItemList *items);
-
-extern double stat_r1_zero(const double *ref1, const double *ref2,
- ReflItemList *items, double *scalep);
-extern double stat_r1_ignore(const double *ref1, const double *ref2,
- ReflItemList *items, double *scalep);
-
-extern double stat_r2(const double *ref1, const double *ref2,
- ReflItemList *items, double *scalep);
-
-extern double stat_r1_i(const double *ref1, const double *ref2,
- ReflItemList *items, double *scalep);
-
-extern double stat_rdiff_zero(const double *ref1, const double *ref2,
- ReflItemList *items, double *scalep);
-extern double stat_rdiff_ignore(const double *ref1, const double *ref2,
- ReflItemList *items, double *scalep);
-extern double stat_rdiff_intensity(const double *ref1, const double *ref2,
- ReflItemList *items, double *scalep);
-
-extern double stat_pearson_i(const double *ref1, const double *ref2,
- ReflItemList *items);
-extern double stat_pearson_f_zero(const double *ref1, const double *ref2,
- ReflItemList *items);
-extern double stat_pearson_f_ignore(const double *ref1, const double *ref2,
- ReflItemList *items);
+#include "reflist.h"
+
+extern double stat_scale_intensity(RefList *list1, RefList *list2);
+
+extern double stat_r1_zero(RefList *list1, RefList *list2, double *scalep);
+extern double stat_r1_ignore(RefList *list1, RefList *list2, double *scalep);
+
+extern double stat_r2(RefList *list1, RefList *list2, double *scalep);
+
+extern double stat_r1_i(RefList *list1, RefList *list2, double *scalep);
+
+extern double stat_rdiff_zero(RefList *list1, RefList *list2, double *scalep);
+extern double stat_rdiff_ignore(RefList *list1, RefList *list2, double *scalep);
+extern double stat_rdiff_intensity(RefList *list1, RefList *list2,
+ double *scalep);
+
+extern double stat_pearson_i(RefList *list1, RefList *list2);
+extern double stat_pearson_f_zero(RefList *list1, RefList *list2);
+extern double stat_pearson_f_ignore(RefList *list1, RefList *list2);
#endif /* STATISTICS_H */
diff --git a/src/stream.c b/src/stream.c
index 545bed43..97c38b77 100644
--- a/src/stream.c
+++ b/src/stream.c
@@ -31,7 +31,8 @@
#define PEAK_LIST_START_MARKER "Peaks from peak search"
#define PEAK_LIST_END_MARKER "End of peak list"
#define REFLECTION_START_MARKER "Reflections measured after indexing"
-#define REFLECTION_END_MARKER "End of reflections"
+/* REFLECTION_END_MARKER is over in reflist-utils.h because it is also
+ * used to terminate a standalone list of reflections */
static void exclusive(const char *a, const char *b)
{
@@ -132,49 +133,6 @@ int count_patterns(FILE *fh)
}
-static int read_reflections(FILE *fh, struct image *image)
-{
- char *rval = NULL;
- int first = 1;
-
- image->reflections = reflist_new();
-
- do {
-
- char line[1024];
- signed int h, k, l;
- float intensity, sigma, res, fs, ss;
- char phs[1024];
- int cts;
- int r;
- Reflection *refl;
-
- rval = fgets(line, 1023, fh);
- if ( rval == NULL ) continue;
- chomp(line);
-
- if ( strcmp(line, REFLECTION_END_MARKER) == 0 ) return 0;
-
- r = sscanf(line, "%i %i %i %f %s %f %f %i %f %f",
- &h, &k, &l, &intensity, phs, &sigma, &res, &cts,
- &fs, &ss);
- if ( (r != 10) && (!first) ) return 1;
-
- first = 0;
- if ( r == 10 ) {
- refl = add_refl(image->reflections, h, k, l);
- set_int(refl, intensity);
- set_detector_pos(refl, fs, ss, 0.0);
- set_esd_intensity(refl, sigma);
- }
-
- } while ( rval != NULL );
-
- /* Got read error of some kind before finding PEAK_LIST_END_MARKER */
- return 1;
-}
-
-
static int read_peaks(FILE *fh, struct image *image)
{
char *rval = NULL;
@@ -409,7 +367,8 @@ int read_chunk(FILE *fh, struct image *image)
}
if ( strcmp(line, REFLECTION_START_MARKER) == 0 ) {
- if ( read_reflections(fh, image) ) {
+ image->reflections = read_reflections_from_file(fh);
+ if ( image->reflections == NULL ) {
ERROR("Failed while reading reflections\n");
return 1;
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-30 09:53:43 +0000