aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/compare_hkl.c69
-rw-r--r--src/diffraction-gpu.c9
-rw-r--r--src/diffraction-gpu.h6
-rw-r--r--src/diffraction.c81
-rw-r--r--src/diffraction.h3
-rw-r--r--src/get_hkl.c94
-rw-r--r--src/indexamajig.c21
-rw-r--r--src/likelihood.c16
-rw-r--r--src/likelihood.h7
-rw-r--r--src/pattern_sim.c15
-rw-r--r--src/process_hkl.c80
-rw-r--r--src/reflections.c83
-rw-r--r--src/reflections.h14
-rw-r--r--src/render_hkl.c10
-rw-r--r--src/sfac.c4
-rw-r--r--src/sfac.h3
-rw-r--r--src/statistics.c197
-rw-r--r--src/statistics.h22
-rw-r--r--src/utils.c42
-rw-r--r--src/utils.h2
20 files changed, 345 insertions, 433 deletions
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index e6caee6f..e69f3753 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -49,13 +49,9 @@ int main(int argc, char *argv[])
char *outfile = NULL;
char *afile = NULL;
char *bfile = NULL;
- signed int h, k, l;
double scale, R2, Rmerge, pearson;
- unsigned int *c1;
- unsigned int *c2;
- int i;
- int nc1, nc2, ncom;
- unsigned int *outcounts;
+ int i, ncom;
+ ReflItemList *i1, *i2, *icommon;
/* Long options */
const struct option longopts[] = {
@@ -94,72 +90,53 @@ int main(int argc, char *argv[])
bfile = strdup(argv[optind]);
cell = load_cell_from_pdb("molecule.pdb");
- c1 = new_list_count();
- ref1 = read_reflections(afile, c1, NULL);
+ ref1 = new_list_intensity();
+ i1 = read_reflections(afile, ref1, NULL, NULL);
if ( ref1 == NULL ) {
ERROR("Couldn't open file '%s'\n", afile);
return 1;
}
- c2 = new_list_count();
- ref2 = read_reflections(bfile, c2, NULL);
+ ref2 = new_list_intensity();
+ i2 = read_reflections(bfile, ref2, NULL, NULL);
if ( ref2 == NULL ) {
ERROR("Couldn't open file '%s'\n", bfile);
return 1;
}
- out = new_list_intensity();
- outcounts = new_list_count();
- /* Knock out the zero beam, in case it's present */
- set_count(c1, 0, 0, 0, 0);
- set_count(c2, 0, 0, 0, 0);
+ /* Find common reflections */
+ icommon = intersection_items(i1, i2);
+ ncom = num_items(icommon);
- /* Divide by number of counts, since we're not interested in them */
- divide_down(ref1, c1);
- divide_down(ref2, c2);
+ /* List for output scale factor map */
+ out = new_list_intensity();
- for ( h=-INDMAX; h<INDMAX; h++ ) {
- for ( k=-INDMAX; k<INDMAX; k++ ) {
- for ( l=-INDMAX; l<INDMAX; l++ ) {
+ for ( i=0; i<ncom; i++ ) {
double i1, i2;
- unsigned int c1s, c2s;
+ struct refl_item *it;
+ signed int h, k, l;
- c1s = lookup_count(c1, h, k, l);
- c2s = lookup_count(c2, h, k, l);
+ it = get_item(icommon, i);
+ h = it->h; k = it->k; l = it->l;
i1 = lookup_intensity(ref1, h, k, l);
i2 = lookup_intensity(ref2, h, k, l);
- if ( c1s && c2s ) {
- if ( (i1 != 0.0) && (i2 != 0.0) ) {
- set_intensity(out, h, k, l,
- (i1/(double)c1s)/i2/(double)c2s);
- set_count(outcounts, h, k, l, 1);
- }
- }
+ set_intensity(out, h, k, l, i1/i2);
}
- }
- }
- nc1 = 0;
- nc2 = 0;
- ncom = 0;
- for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
- nc1 += c1[i];
- nc2 += c2[i];
- ncom += c1[i] && c2[i];
- }
- STATUS("%i,%i reflections: %i in common\n", nc1, nc2, ncom);
- R2 = stat_r2(ref1, c1, ref2, c2, &scale);
+ STATUS("%i,%i reflections: %i in common\n",
+ num_items(i1), num_items(i2), ncom);
+ R2 = stat_r2(ref1, ref2, icommon, &scale);
STATUS("R2 = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale);
- Rmerge = stat_rmerge(ref1, c1, ref2, c2, &scale);
+ Rmerge = stat_rmerge(ref1, ref2, icommon, &scale);
STATUS("Rmerge = %5.4f %% (scale=%5.2e)\n", Rmerge*100.0, scale);
- pearson = stat_pearson(ref1, c1, ref2, c2);
+ pearson = stat_pearson(ref1, ref2, icommon);
STATUS("Pearson r = %5.4f\n", pearson);
if ( outfile != NULL ) {
- write_reflections(outfile, outcounts, out, NULL, 0, cell, 1);
+ write_reflections(outfile, icommon, out, NULL, NULL, cell);
}
return 0;
diff --git a/src/diffraction-gpu.c b/src/diffraction-gpu.c
index 86594a47..7a5bef85 100644
--- a/src/diffraction-gpu.c
+++ b/src/diffraction-gpu.c
@@ -342,8 +342,7 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
/* Setup the OpenCL stuff, create buffers, load the structure factor table */
struct gpu_context *setup_gpu(int no_sfac, struct image *image,
- const double *intensities,
- const unsigned int *counts)
+ const double *intensities)
{
struct gpu_context *gctx;
cl_uint nplat;
@@ -404,11 +403,7 @@ struct gpu_context *setup_gpu(int no_sfac, struct image *image,
intensities_ptr = malloc(intensities_size);
if ( intensities != NULL ) {
for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
- if ( counts[i] == 1 ) {
- intensities_ptr[i] = intensities[i];
- } else {
- intensities_ptr[i] = 1.0e20;
- }
+ intensities_ptr[i] = intensities[i];
}
} else {
for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
diff --git a/src/diffraction-gpu.h b/src/diffraction-gpu.h
index 003c4c3a..2438e8f1 100644
--- a/src/diffraction-gpu.h
+++ b/src/diffraction-gpu.h
@@ -26,8 +26,7 @@ struct gpu_context;
extern void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
int na, int nb, int nc, UnitCell *ucell);
extern struct gpu_context *setup_gpu(int no_sfac, struct image *image,
- const double *intensities,
- const unsigned int *counts);
+ const double *intensities);
extern void cleanup_gpu(struct gpu_context *gctx);
#else
@@ -40,8 +39,7 @@ static void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
}
static struct gpu_context *setup_gpu(int no_sfac, struct image *image,
- const double *intensities,
- const unsigned int *counts)
+ const double *intensities)
{
return NULL;
}
diff --git a/src/diffraction.c b/src/diffraction.c
index fb512e62..efc56231 100644
--- a/src/diffraction.c
+++ b/src/diffraction.c
@@ -71,13 +71,11 @@ static double lattice_factor(struct rvec q, double ax, double ay, double az,
static double interpolate_linear(const double *ref,
- const unsigned int *counts,
float hd, signed int k, signed int l)
{
signed int h;
double val1, val2;
float f;
- unsigned int c1, c2;
h = (signed int)hd;
if ( hd < 0.0 ) h -= 1;
@@ -86,30 +84,15 @@ static double interpolate_linear(const double *ref,
val1 = lookup_intensity(ref, h, k, l);
val2 = lookup_intensity(ref, h+1, k, l);
- c1 = lookup_count(counts, h, k, l);
- c2 = lookup_count(counts, h+1, k, l);
- if ( c1 == 0 ) {
- ERROR("Needed intensity for %i %i %i, but don't have it.\n",
- h, k, l);
- return 1.0e20;
- }
-
- if ( c2 == 0 ) {
- ERROR("Needed intensity for %i %i %i, but don't have it.\n",
- h+1, k, l);
- return 1.0e20;
- }
-
- val1 = val1 / (double)c1;
- val2 = val2 / (double)c2;
+ val1 = val1;
+ val2 = val2;
return (1.0-f)*val1 + f*val2;
}
static double interpolate_bilinear(const double *ref,
- const unsigned int *counts,
float hd, float kd, signed int l)
{
signed int k;
@@ -121,15 +104,14 @@ static double interpolate_bilinear(const double *ref,
f = kd - (float)k;
assert(f >= 0.0);
- val1 = interpolate_linear(ref, counts, hd, k, l);
- val2 = interpolate_linear(ref, counts, hd, k+1, l);
+ val1 = interpolate_linear(ref, hd, k, l);
+ val2 = interpolate_linear(ref, hd, k+1, l);
return (1.0-f)*val1 + f*val2;
}
static double interpolate_intensity(const double *ref,
- const unsigned int *counts,
float hd, float kd, float ld)
{
signed int l;
@@ -141,15 +123,14 @@ static double interpolate_intensity(const double *ref,
f = ld - (float)l;
assert(f >= 0.0);
- val1 = interpolate_bilinear(ref, counts, hd, kd, l);
- val2 = interpolate_bilinear(ref, counts, hd, kd, l+1);
+ val1 = interpolate_bilinear(ref, hd, kd, l);
+ val2 = interpolate_bilinear(ref, hd, kd, l+1);
return (1.0-f)*val1 + f*val2;
}
static double complex interpolate_phased_linear(const double *ref,
- const unsigned int *counts,
const double *phases,
float hd,
signed int k, signed int l)
@@ -157,7 +138,6 @@ static double complex interpolate_phased_linear(const double *ref,
signed int h;
double val1, val2;
float f;
- unsigned int c1, c2;
double ph1, ph2;
double re1, re2, im1, im2;
double re, im;
@@ -169,25 +149,11 @@ static double complex interpolate_phased_linear(const double *ref,
val1 = lookup_intensity(ref, h, k, l);
val2 = lookup_intensity(ref, h+1, k, l);
- c1 = lookup_count(counts, h, k, l);
- c2 = lookup_count(counts, h+1, k, l);
ph1 = lookup_phase(phases, h, k, l);
ph2 = lookup_phase(phases, h+1, k, l);
- if ( c1 == 0 ) {
- ERROR("Needed intensity for %i %i %i, but don't have it.\n",
- h, k, l);
- return 1.0e20;
- }
-
- if ( c2 == 0 ) {
- ERROR("Needed intensity for %i %i %i, but don't have it.\n",
- h+1, k, l);
- return 1.0e20;
- }
-
- val1 = val1 / (double)c1;
- val2 = val2 / (double)c2;
+ val1 = val1;
+ val2 = val2;
/* Calculate real and imaginary parts */
re1 = val1 * cos(ph1);
@@ -203,7 +169,6 @@ static double complex interpolate_phased_linear(const double *ref,
static double complex interpolate_phased_bilinear(const double *ref,
- const unsigned int *counts,
const double *phases,
float hd, float kd,
signed int l)
@@ -217,15 +182,14 @@ static double complex interpolate_phased_bilinear(const double *ref,
f = kd - (float)k;
assert(f >= 0.0);
- val1 = interpolate_phased_linear(ref, counts, phases, hd, k, l);
- val2 = interpolate_phased_linear(ref, counts, phases, hd, k+1, l);
+ val1 = interpolate_phased_linear(ref, phases, hd, k, l);
+ val2 = interpolate_phased_linear(ref, phases, hd, k+1, l);
return (1.0-f)*val1 + f*val2;
}
static double interpolate_phased_intensity(const double *ref,
- const unsigned int *counts,
const double *phases,
float hd, float kd, float ld)
{
@@ -238,16 +202,15 @@ static double interpolate_phased_intensity(const double *ref,
f = ld - (float)l;
assert(f >= 0.0);
- val1 = interpolate_phased_bilinear(ref, counts, phases, hd, kd, l);
- val2 = interpolate_phased_bilinear(ref, counts, phases, hd, kd, l+1);
+ val1 = interpolate_phased_bilinear(ref, phases, hd, kd, l);
+ val2 = interpolate_phased_bilinear(ref, phases, hd, kd, l+1);
return cabs((1.0-f)*val1 + f*val2);
}
/* Look up the structure factor for the nearest Bragg condition */
-static double molecule_factor(const double *intensities,
- const unsigned int *counts, const double *phases,
+static double molecule_factor(const double *intensities,const double *phases,
struct rvec q,
double ax, double ay, double az,
double bx, double by, double bz,
@@ -267,18 +230,13 @@ static double molecule_factor(const double *intensities,
h = (signed int)rint(hd);
k = (signed int)rint(kd);
l = (signed int)rint(ld);
- if ( lookup_count(counts, h, k, l) == 0 ) {
- ERROR("Needed intensity for %i %i %i, but don't have it.\n",
- h, k, l);
- return 1.0e20;
- }
r = lookup_intensity(intensities, h, k, l);
break;
case GRADIENT_INTERPOLATE :
- r = interpolate_intensity(intensities, counts, hd, kd, ld);
+ r = interpolate_intensity(intensities, hd, kd, ld);
break;
case GRADIENT_PHASED :
- r = interpolate_phased_intensity(intensities, counts, phases,
+ r = interpolate_phased_intensity(intensities, phases,
hd, kd, ld);
break;
default:
@@ -291,7 +249,7 @@ static double molecule_factor(const double *intensities,
double water_diffraction(struct rvec q, double en,
- double beam_r, double water_r)
+ double beam_r, double water_r)
{
double complex fH, fO;
double s, modq;
@@ -390,9 +348,8 @@ double get_tt(struct image *image, unsigned int xs, unsigned int ys)
void get_diffraction(struct image *image, int na, int nb, int nc,
- const double *intensities, const unsigned int *counts,
- const double *phases, UnitCell *cell, int do_water,
- GradientMethod m)
+ const double *intensities, const double *phases,
+ UnitCell *cell, int do_water, GradientMethod m)
{
unsigned int xs, ys;
double ax, ay, az;
@@ -447,7 +404,7 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
I_molecule = 1.0e10;
} else {
I_molecule = molecule_factor(intensities,
- counts, phases, q,
+ phases, q,
ax,ay,az,
bx,by,bz,cx,cy,cz,
m);
diff --git a/src/diffraction.h b/src/diffraction.h
index c7afaa3e..ffe3a19d 100644
--- a/src/diffraction.h
+++ b/src/diffraction.h
@@ -26,8 +26,7 @@ typedef enum {
} GradientMethod;
extern void get_diffraction(struct image *image, int na, int nb, int nc,
- const double *intensities,
- const unsigned int *counts, const double *phases,
+ const double *intensities,const double *phases,
UnitCell *cell, int do_water, GradientMethod m);
extern struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
unsigned int sampling, float *ttp, float k);
diff --git a/src/get_hkl.c b/src/get_hkl.c
index c98508ec..78d79f76 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -38,8 +38,6 @@ static void show_help(const char *s)
" --poisson Simulate Poisson samples.\n"
" --twin Generate twinned data.\n"
" -o, --output=<filename> Output filename (default: stdout).\n"
-" --zone-axis Generate hk0 intensities only (and add\n"
-" Synth2D-style header.\n"
" -i, --intensities=<file> Read intensities from file instead of\n"
" calculating them from scratch. You might use\n"
" this if you need to apply noise or twinning.\n"
@@ -48,37 +46,6 @@ static void show_help(const char *s)
}
-static int template_reflections(const char *filename, unsigned int *counts)
-{
- char *rval;
- FILE *fh;
-
- fh = fopen(filename, "r");
- if ( fh == NULL ) {
- return 1;
- }
-
- do {
-
- char line[1024];
- int r;
- signed int h, k, l;
-
- rval = fgets(line, 1023, fh);
-
- r = sscanf(line, "%i %i %i", &h, &k, &l);
- if ( r != 3 ) continue;
-
- set_count(counts, h, k, l, 1);
-
- } while ( rval != NULL );
-
- fclose(fh);
-
- return 0;
-}
-
-
/* Apply Poisson noise to all reflections */
static void noisify_reflections(double *ref)
{
@@ -111,13 +78,12 @@ int main(int argc, char *argv[])
char *template = NULL;
int config_noisify = 0;
int config_twin = 0;
- int config_za = 0;
char *output = NULL;
- unsigned int *counts;
- unsigned int *cts;
char *input = NULL;
signed int h, k, l;
char *filename = NULL;
+ ReflItemList *input_items;
+ ReflItemList *write_items;
/* Long options */
const struct option longopts[] = {
@@ -126,7 +92,6 @@ int main(int argc, char *argv[])
{"poisson", 0, &config_noisify, 1},
{"output", 1, NULL, 'o'},
{"twin", 0, &config_twin, 1},
- {"zone-axis", 0, &config_za, 1},
{"intensities", 1, NULL, 'i'},
{"pdb", 1, NULL, 'p'},
{0, 0, NULL, 0}
@@ -170,41 +135,18 @@ int main(int argc, char *argv[])
}
mol = load_molecule(filename);
- cts = new_list_count();
- phases = new_list_intensity(); /* "intensity" type used for phases */
+ phases = new_list_phase();
if ( input == NULL ) {
+ input_items = new_items();
ideal_ref = get_reflections(mol, eV_to_J(1790.0), 1/(0.05e-9),
- cts, phases);
+ phases, input_items);
} else {
- ideal_ref = read_reflections(input, cts, phases);
+ ideal_ref = new_list_intensity();
+ phases = new_list_phase();
+ input_items = read_reflections(input, ideal_ref, phases, NULL);
free(input);
}
- counts = new_list_count();
-
- if ( template != NULL ) {
-
- if ( template_reflections(template, counts) != 0 ) {
- ERROR("Failed to template reflections.\n");
- return 1;
- }
-
- } else {
-
- /* No template? Then only mark reflections which were
- * calculated. */
- for ( h=-INDMAX; h<=INDMAX; h++ ) {
- for ( k=-INDMAX; k<=INDMAX; k++ ) {
- for ( l=-INDMAX; l<=INDMAX; l++ ) {
- unsigned int c;
- c = lookup_count(cts, h, k, l);
- set_count(counts, h, k, l, c);
- }
- }
- }
-
- }
-
if ( config_noisify ) noisify_reflections(ideal_ref);
if ( config_twin ) {
@@ -239,8 +181,24 @@ int main(int argc, char *argv[])
}
- write_reflections(output, counts, ideal_ref, phases,
- config_za, mol->cell, 1);
+ 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);
+ write_items = intersection_items(input_items, template_items);
+ delete_items(template_items);
+ } else {
+ /* Write out all reflections */
+ write_items = new_items();
+ union_items(write_items, input_items);
+ }
+
+ write_reflections(output, write_items, ideal_ref, phases, NULL,
+ mol->cell);
+
+ delete_items(input_items);
+ delete_items(write_items);
return 0;
}
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 3a3dfb56..5df08ad3 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -65,7 +65,6 @@ struct process_args
struct detector *det;
IndexingMethod indm;
const double *intensities;
- const unsigned int *counts;
struct gpu_context *gctx;
/* Thread control and output */
@@ -225,19 +224,18 @@ static struct image *get_simage(struct image *template, int alternate)
static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
- const double *intensities,
- const unsigned int *counts, UnitCell *cell)
+ const double *intensities, UnitCell *cell)
{
/* Set up GPU if necessary */
if ( (gctx != NULL) && (*gctx == NULL) ) {
- *gctx = setup_gpu(0, simage, intensities, counts);
+ *gctx = setup_gpu(0, simage, intensities);
}
if ( (gctx != NULL) && (*gctx != NULL) ) {
get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
} else {
get_diffraction(simage, 24, 24, 40,
- intensities, counts, NULL, cell, 0,
+ intensities, NULL, cell, 0,
GRADIENT_MOSAIC);
}
record_image(simage, 0);
@@ -270,7 +268,6 @@ static struct process_result process_image(struct process_args *pargs)
int config_polar = pargs->config_polar;
IndexingMethod indm = pargs->indm;
const double *intensities = pargs->intensities;
- const unsigned int *counts = pargs->counts;
struct gpu_context *gctx = pargs->gctx;
image.features = NULL;
@@ -368,11 +365,11 @@ static struct process_result process_image(struct process_args *pargs)
if ( config_gpu ) {
pthread_mutex_lock(pargs->gpu_mutex);
simulate_and_write(simage, &gctx, intensities,
- counts, image.indexed_cell);
+ image.indexed_cell);
pthread_mutex_unlock(pargs->gpu_mutex);
} else {
simulate_and_write(simage, NULL, intensities,
- counts, image.indexed_cell);
+ image.indexed_cell);
}
}
@@ -469,7 +466,6 @@ int main(int argc, char *argv[])
UnitCell *cell;
double *intensities = NULL;
char *intfile = NULL;
- unsigned int *counts;
char *pdb = NULL;
char *prefix = NULL;
int nthreads = 1;
@@ -568,11 +564,11 @@ int main(int argc, char *argv[])
free(filename);
if ( intfile != NULL ) {
- counts = new_list_count();
- intensities = read_reflections(intfile, counts, NULL);
+ ReflItemList *items;
+ items = read_reflections(intfile, intensities, NULL, NULL);
+ delete_items(items);
} else {
intensities = NULL;
- counts = NULL;
}
if ( pdb == NULL ) {
@@ -671,7 +667,6 @@ int main(int argc, char *argv[])
pargs->det = det;
pargs->indm = indm;
pargs->intensities = intensities;
- pargs->counts = counts;
pargs->gctx = gctx;
pargs->id = i;
pthread_mutex_lock(&pargs->control_mutex);
diff --git a/src/likelihood.c b/src/likelihood.c
index b6a37994..cb5595f8 100644
--- a/src/likelihood.c
+++ b/src/likelihood.c
@@ -18,17 +18,17 @@
#include "utils.h"
-void scale_intensities(const double *model, double *new_pattern,
- const unsigned int *model_counts,
- ReflItemList *items, double f0, int f0_valid)
+void scale_intensities(const double *model, ReflItemList *model_items,
+ double *new_pattern, ReflItemList *new_items,
+ double f0, int f0_valid)
{
double s;
unsigned int i;
- unsigned int *new_counts;
+ ReflItemList *items;
- new_counts = items_to_counts(items);
-
- s = stat_scale_intensity(model, model_counts, new_pattern, new_counts);
+ items = intersection_items(model_items, new_items);
+ s = stat_scale_intensity(model, new_pattern, items);
+ delete_items(items);
if ( f0_valid ) printf("%f %f\n", s, f0);
/* NaN -> abort */
@@ -38,6 +38,4 @@ void scale_intensities(const double *model, double *new_pattern,
for ( i=0; i<LIST_SIZE; i++ ) {
new_pattern[i] *= s;
}
-
- free(new_counts);
}
diff --git a/src/likelihood.h b/src/likelihood.h
index 08956fe3..70aef098 100644
--- a/src/likelihood.h
+++ b/src/likelihood.h
@@ -20,10 +20,9 @@
#include "utils.h"
-extern void scale_intensities(const double *model, double *new_pattern,
- const unsigned int *model_counts,
- ReflItemList *items, double f0,
- int f0_valid);
+extern void scale_intensities(const double *model, ReflItemList *model_items,
+ double *new_pattern, ReflItemList *new_items,
+ double f0, int f0_valid);
#endif /* LIKELIHOOD_H */
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index 5d59d71b..616eec2c 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -195,7 +195,6 @@ int main(int argc, char *argv[])
int n_images = 1; /* Generate one image by default */
int done = 0;
UnitCell *cell;
- unsigned int *counts;
/* Long options */
const struct option longopts[] = {
@@ -323,17 +322,19 @@ int main(int argc, char *argv[])
STATUS("reflection intensities (with --intensities).\n");
STATUS("I'll simulate a flat intensity distribution.\n");
intensities = NULL;
- counts = NULL;
phases = NULL;
} else {
- counts = new_list_count();
+ ReflItemList *items;
if ( grad == GRADIENT_PHASED ) {
phases = new_list_phase();
} else {
phases = NULL;
}
- intensities = read_reflections(intfile, counts, phases);
+ intensities = new_list_intensity();
+ phases = new_list_phase();
+ items = read_reflections(intfile, intensities, phases, NULL);
free(intfile);
+ delete_items(items);
}
/* Define image parameters */
@@ -401,12 +402,12 @@ int main(int argc, char *argv[])
if ( config_gpu ) {
if ( gctx == NULL ) {
gctx = setup_gpu(config_nosfac, &image,
- intensities, counts);
+ intensities);
}
get_diffraction_gpu(gctx, &image, na, nb, nc, cell);
} else {
- get_diffraction(&image, na, nb, nc, intensities, counts,
- phases, cell, !config_nowater, grad);
+ get_diffraction(&image, na, nb, nc, intensities, phases,
+ cell, !config_nowater, grad);
}
if ( image.data == NULL ) {
ERROR("Diffraction calculation failed.\n");
diff --git a/src/process_hkl.c b/src/process_hkl.c
index b5b6df9a..986477fb 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -89,6 +89,8 @@ static ReflItemList *get_twin_possibilities(const char *holo, const char *mero)
twins = get_twins(test_items, holo, mero);
delete_items(test_items);
+ /* Idiot check. Wouldn't be necessary if I could prove that the above
+ * set of arbitrarily chosen reflections were always general. */
if ( num_items(twins) != np ) {
ERROR("Whoops! Couldn't find all the twinning possiblities.\n");
abort();
@@ -98,10 +100,9 @@ static ReflItemList *get_twin_possibilities(const char *holo, const char *mero)
}
-static int resolve_twin(ReflItemList *items, ReflItemList *twins,
- const double *patt, const double *model,
- const unsigned int *model_counts,
- const char *holo, const char *mero)
+static int resolve_twin(const double *model, ReflItemList *observed,
+ const double *patt, ReflItemList *items,
+ ReflItemList *twins, const char *holo, const char *mero)
{
int n, i;
double best_fom = 0.0;
@@ -116,6 +117,7 @@ static int resolve_twin(ReflItemList *items, ReflItemList *twins,
double *trial_ints = new_list_intensity();
unsigned int *trial_counts = new_list_count();
double fom;
+ ReflItemList *intersection;
op = get_item(twins, i)->op;
@@ -134,8 +136,9 @@ static int resolve_twin(ReflItemList *items, ReflItemList *twins,
}
- fom = stat_pearson(trial_ints, trial_counts,
- model, model_counts);
+ intersection = intersection_items(observed, items);
+ fom = stat_pearson(trial_ints, model, intersection);
+ delete_items(intersection);
free(trial_ints);
free(trial_counts);
@@ -153,9 +156,9 @@ static int resolve_twin(ReflItemList *items, ReflItemList *twins,
}
-static void merge_pattern(double *model, const double *new,
- unsigned int *model_counts,
- ReflItemList *items, int mo, int sum,
+static void merge_pattern(double *model, ReflItemList *observed,
+ const double *new, ReflItemList *items,
+ unsigned int *model_counts, int mo,
ReflItemList *twins,
const char *holo, const char *mero)
{
@@ -163,8 +166,8 @@ static void merge_pattern(double *model, const double *new,
int twin;
if ( twins != NULL ) {
- twin = resolve_twin(items, twins, new, model,
- model_counts, holo, mero);
+ twin = resolve_twin(model, observed, new, items,
+ twins, holo, mero);
} else {
twin = 0;
}
@@ -187,25 +190,23 @@ static void merge_pattern(double *model, const double *new,
intensity = lookup_intensity(new, h, k, l);
+ /* User asked for max only? */
if ( !mo ) {
integrate_intensity(model, h, k, l, intensity);
- if ( sum ) {
- set_count(model_counts, h, k, l, 1);
- } else {
- integrate_count(model_counts, h, k, l, 1);
- }
} else {
if ( intensity > lookup_intensity(model, h, k, l) ) {
set_intensity(model, h, k, l, intensity);
}
- set_count(model_counts, h, k, l, 1);
}
+ integrate_count(model_counts, h, k, l, 1);
+
}
}
-static void merge_all(FILE *fh, double *model, unsigned int *model_counts,
+static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved,
+ unsigned int **pcounts,
int config_maxonly, int config_scale, int config_sum,
int config_stopafter,
ReflItemList *twins, const char *holo, const char *mero,
@@ -218,6 +219,10 @@ static void merge_all(FILE *fh, double *model, unsigned int *model_counts,
int n_patterns = 0;
double *new_pattern = new_list_intensity();
ReflItemList *items = new_items();
+ ReflItemList *observed = new_items();
+ double *model = new_list_intensity();
+ unsigned int *counts = new_list_count();
+ int i;
do {
@@ -244,19 +249,20 @@ static void merge_all(FILE *fh, double *model, unsigned int *model_counts,
/* Scale if requested */
if ( config_scale ) {
- scale_intensities(model, new_pattern,
- model_counts, items, f0,
- f0_valid);
+ scale_intensities(model, observed,
+ new_pattern, items,
+ f0, f0_valid);
}
/* Start of second or later pattern */
- merge_pattern(model, new_pattern, model_counts,
- items, config_maxonly, config_sum, twins,
- holo, mero);
+ merge_pattern(model, observed, new_pattern, items,
+ counts, config_maxonly,
+ twins, holo, mero);
if ( n_patterns == config_stopafter ) break;
n_patterns++;
+ union_items(observed, items);
clear_items(items);
progress_bar(n_patterns, n_total_patterns, "Merging");
@@ -292,6 +298,19 @@ static void merge_all(FILE *fh, double *model, unsigned int *model_counts,
delete_items(items);
free(new_pattern);
+ if ( !config_sum ) {
+ for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
+ if ( counts[i] > 0 ) {
+ model[i] /= (double)counts[i];
+ counts[i] = 1;
+ }
+ }
+ }
+
+ *pmodel = model;
+ *pcounts = counts;
+ *pobserved = observed;
+
STATUS("%i patterns had no f0 valid value.\n", n_nof0);
}
@@ -322,7 +341,7 @@ int main(int argc, char *argv[])
char *output = NULL;
FILE *fh;
double *model;
- unsigned int *model_counts;
+ unsigned int *counts;
UnitCell *cell;
int config_maxonly = 0;
int config_stopafter = 0;
@@ -332,6 +351,7 @@ int main(int argc, char *argv[])
char *sym = NULL;
char *pdb = NULL;
ReflItemList *twins;
+ ReflItemList *observed;
int i;
/* Long options */
@@ -400,9 +420,6 @@ int main(int argc, char *argv[])
sym = strdup("1");
}
- model = new_list_intensity();
- model_counts = new_list_count();
-
cell = load_cell_from_pdb(pdb);
free(pdb);
@@ -443,7 +460,7 @@ int main(int argc, char *argv[])
STATUS("There are %i patterns to process\n", n_total_patterns);
rewind(fh);
- merge_all(fh, model, model_counts,
+ merge_all(fh, &model, &observed, &counts,
config_maxonly, config_scale, config_sum, config_stopafter,
twins, holo, sym, n_total_patterns);
rewind(fh);
@@ -451,13 +468,12 @@ int main(int argc, char *argv[])
fclose(fh);
if ( output != NULL ) {
- write_reflections(output, model_counts, model, NULL,
- 0, cell, 1);
+ write_reflections(output, observed, model, NULL, counts, cell);
}
free(sym);
free(model);
- free(model_counts);
+ free(counts);
free(output);
free(cell);
diff --git a/src/reflections.c b/src/reflections.c
index bd9b446d..0ab7a91e 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -21,12 +21,12 @@
#include "reflections.h"
-void write_reflections(const char *filename, unsigned int *counts,
- double *ref, double *phases, int zone_axis,
- UnitCell *cell, unsigned int min_counts)
+void write_reflections(const char *filename, ReflItemList *items,
+ double *intensities, double *phases,
+ unsigned int *counts, UnitCell *cell)
{
FILE *fh;
- signed int h, k, l;
+ int i;
if ( filename == NULL ) {
fh = stdout;
@@ -40,23 +40,16 @@ void write_reflections(const char *filename, unsigned int *counts,
}
/* Write spacings and angle if zone axis pattern */
- if ( zone_axis ) {
- double a, b, c, alpha, beta, gamma;
- cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma);
- fprintf(fh, "a %5.3f nm\n",
- (0.5/resolution(cell, 1, 0, 0))*1e9);
- fprintf(fh, "b %5.3f nm\n",
- (0.5/resolution(cell, 0, 1, 0))*1e9);
- fprintf(fh, "angle %5.3f deg\n", rad2deg(alpha));
- fprintf(fh, "scale 10\n");
- } else {
- fprintf(fh, " h k l I phase sigma(I) "
- " 1/d(nm^-1) counts\n");
- }
+ 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;
- for ( h=-INDMAX; h<INDMAX; h++ ) {
- for ( k=-INDMAX; k<INDMAX; k++ ) {
- for ( l=-INDMAX; l<INDMAX; l++ ) {
+ it = get_item(items, i);
+ h = it->h; k = it->k; l = it->l;
int N;
double intensity, s;
@@ -64,14 +57,12 @@ void write_reflections(const char *filename, unsigned int *counts,
if ( counts ) {
N = lookup_count(counts, h, k, l);
- if ( N < min_counts ) continue;
} else {
N = 1;
}
- if ( zone_axis && (l != 0) ) continue;
- /* Divide measured intensity by the number of counts */
- intensity = lookup_intensity(ref, h, k, l) / N;
+ intensity = lookup_intensity(intensities, h, k, l);
+
if ( phases != NULL ) {
double p;
p = lookup_phase(phases, h, k, l);
@@ -91,18 +82,25 @@ void write_reflections(const char *filename, unsigned int *counts,
h, k, l, intensity, ph, 0.0, s/1.0e9, N);
}
- }
- }
fclose(fh);
}
-double *read_reflections(const char *filename, unsigned int *counts,
- double *phases)
+/* 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 *ref;
FILE *fh;
char *rval;
+ ReflItemList *items;
fh = fopen(filename, "r");
if ( fh == NULL ) {
@@ -110,7 +108,7 @@ double *read_reflections(const char *filename, unsigned int *counts,
return NULL;
}
- ref = new_list_intensity();
+ items = new_items();
do {
@@ -147,24 +145,24 @@ double *read_reflections(const char *filename, unsigned int *counts,
continue;
}
- set_intensity(ref, h, k, l, intensity);
+ 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);
- /* In this case, the intensity must be multiplied up
- * because other parts of the program will try to
- * divide it down. */
- set_intensity(ref, h, k, l, intensity*(double)cts);
}
} while ( rval != NULL );
fclose(fh);
- return ref;
+ return items;
}
@@ -188,16 +186,3 @@ double *ideal_intensities(double complex *sfac)
return ref;
}
-
-
-void divide_down(double *intensities, unsigned int *counts)
-{
- int i;
-
- for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
- if ( counts[i] > 0 ) {
- intensities[i] /= (double)counts[i];
- counts[i] = 1;
- }
- }
-}
diff --git a/src/reflections.h b/src/reflections.h
index 84946f1a..de5feb72 100644
--- a/src/reflections.h
+++ b/src/reflections.h
@@ -18,18 +18,18 @@
#include "cell.h"
+#include "utils.h"
-extern void write_reflections(const char *filename, unsigned int *counts,
- double *ref, double *phases, int zone_axis,
- UnitCell *cell, unsigned int min_counts);
+extern void write_reflections(const char *filename, ReflItemList *items,
+ double *intensities, double *phases,
+ unsigned int *counts, UnitCell *cell);
-extern double *read_reflections(const char *filename, unsigned int *counts,
- double *phases);
+extern ReflItemList *read_reflections(const char *filename,
+ double *intensities, double *phases,
+ unsigned int *counts);
extern double *ideal_intensities(double complex *sfac);
-extern void divide_down(double *intensities, unsigned int *counts);
-
#endif /* REFLECTIONS_H */
diff --git a/src/render_hkl.c b/src/render_hkl.c
index 60153652..829ceb4e 100644
--- a/src/render_hkl.c
+++ b/src/render_hkl.c
@@ -128,10 +128,10 @@ static void render_za(UnitCell *cell, double *ref, unsigned int *c,
switch ( wght ) {
case WGHT_I :
- val = lookup_intensity(ref, h, k, 0) / (float)ct;
+ val = lookup_intensity(ref, h, k, 0);
break;
case WGHT_SQRTI :
- val = lookup_intensity(ref, h, k, 0) / (float)ct;
+ val = lookup_intensity(ref, h, k, 0);
val = (val>0.0) ? sqrt(val) : 0.0;
break;
case WGHT_COUNTS :
@@ -290,7 +290,6 @@ int main(int argc, char *argv[])
UnitCell *cell;
char *infile;
double *ref;
- unsigned int *cts;
int config_povray = 0;
int config_zoneaxis = 0;
int config_sqrt = 0;
@@ -303,6 +302,7 @@ int main(int argc, char *argv[])
int wght;
int colscale;
char *cscale = NULL;
+ unsigned int *cts;
/* Long options */
const struct option longopts[] = {
@@ -417,8 +417,10 @@ 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();
- ref = read_reflections(infile, cts, NULL);
+ ReflItemList *items = read_reflections(infile, ref, NULL, cts);
+ delete_items(items);
if ( ref == NULL ) {
ERROR("Couldn't open file '%s'\n", infile);
return 1;
diff --git a/src/sfac.c b/src/sfac.c
index a8b94051..b5ae7ace 100644
--- a/src/sfac.c
+++ b/src/sfac.c
@@ -491,7 +491,7 @@ void free_molecule(struct molecule *mol)
double *get_reflections(struct molecule *mol, double en, double res,
- unsigned int *counts, double *phases)
+ double *phases, ReflItemList *items)
{
double *reflections;
double asx, asy, asz;
@@ -556,7 +556,7 @@ double *get_reflections(struct molecule *mol, double en, double res,
set_intensity(reflections, h, k, l, pow(cabs(F), 2.0));
if ( phases != NULL ) set_phase(phases, h, k, l, carg(F));
- set_count(counts, h, k, l, 1);
+ if ( items != NULL ) add_item(items, h, k, l);
}
progress_bar((k+INDMAX)+IDIM*(h+INDMAX), IDIM*IDIM-1,
diff --git a/src/sfac.h b/src/sfac.h
index ea30c83e..115dcb5f 100644
--- a/src/sfac.h
+++ b/src/sfac.h
@@ -19,6 +19,7 @@
#include <complex.h>
#include "cell.h"
+#include "utils.h"
#define MAX_ATOMS (128*1024)
@@ -61,7 +62,7 @@ extern double complex get_sfac(const char *n, double s, double en);
extern struct molecule *load_molecule(const char *filename);
extern void free_molecule(struct molecule *mol);
extern double *get_reflections(struct molecule *mol, double en, double res,
- unsigned int *counts, double *phases);
+ double *phases, ReflItemList *items);
#endif /* SFAC_H */
diff --git a/src/statistics.c b/src/statistics.c
index 889ade51..1acbb718 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -26,10 +26,9 @@
struct r_params {
const double *ref1;
- const unsigned int *c1;
const double *ref2;
- const unsigned int *c2;
- int fom;
+ ReflItemList *items; /* Which reflections to use */
+ int fom; /* Which FoM to use (see the enum just below) */
};
enum {
@@ -38,23 +37,31 @@ enum {
};
-double stat_scale_intensity(const double *ref1, const unsigned int *c1,
- const double *ref2, const unsigned int *c2)
+/* Return the least squares optimal scaling factor when comparing intensities.
+ * ref1,ref2 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 top = 0.0;
double bot = 0.0;
int i;
- /* Start from i=1 -> skip central beam */
- for ( i=1; i<LIST_SIZE; i++ ) {
+ for ( i=0; i<num_items(items); i++ ) {
+
+ double i1, i2;
+ struct refl_item *it;
+ signed int h, k, l;
+
+ it = get_item(items, i);
+ h = it->h; k = it->k; l = it->l;
+
+ i1 = lookup_intensity(ref1, h, k, l);
+ i2 = lookup_intensity(ref2, h, k, l);
- if ( c1[i] && c2[i] ) {
- double i1, i2;
- i1 = ref1[i] / (double)c1[i];
- i2 = ref2[i] / (double)c2[i];
- top += i1 * i2;
- bot += i2 * i2;
- } /* else reflection not common so don't worry about it */
+ top += i1 * i2;
+ bot += i2 * i2;
}
@@ -62,29 +69,36 @@ double stat_scale_intensity(const double *ref1, const unsigned int *c1,
}
-double stat_scale_sqrti(const double *ref1, const unsigned int *c1,
- const double *ref2, const unsigned int *c2)
+/* 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,
+ * 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)
{
double top = 0.0;
double bot = 0.0;
int i;
- /* Start from i=1 -> skip central beam */
- for ( i=1; i<LIST_SIZE; i++ ) {
-
- if ( c1[i] && c2[i] ) {
+ for ( i=0; i<num_items(items); i++ ) {
- double f1, f2;
+ double i1, i2, f1, f2;
+ struct refl_item *it;
+ signed int h, k, l;
- if ( (ref1[i]<0.0) || (ref2[i]<0.0) ) continue;
+ it = get_item(items, i);
+ h = it->h; k = it->k; l = it->l;
- f1 = sqrt(ref1[i]) / (double)c1[i];
- f2 = sqrt(ref2[i]) / (double)c2[i];
-
- top += f1 * f2;
- bot += f2 * f2;
+ 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;
- } /* else reflection not common so don't worry about it */
+ top += i1 * i2;
+ bot += i2 * i2;
}
@@ -92,27 +106,27 @@ double stat_scale_sqrti(const double *ref1, const unsigned int *c1,
}
-static double internal_r2(const double *ref1, const unsigned int *c1,
- const double *ref2, const unsigned int *c2,
- double scale)
+static double internal_r2(const double *ref1, const double *ref2,
+ ReflItemList *items, double scale)
{
double top = 0.0;
double bot = 0.0;
int i;
- /* Start from i=1 -> skip central beam */
- for ( i=1; i<LIST_SIZE; i++ ) {
+ for ( i=0; i<num_items(items); i++ ) {
- if ( c1[i] && c2[i] ) {
+ double i1, i2;
+ struct refl_item *it;
+ signed int h, k, l;
- double i1, i2;
- i1 = ref1[i] / (scale*(double)c1[i]);
- i2 = ref2[i] / (double)c2[i];
+ it = get_item(items, i);
+ h = it->h; k = it->k; l = it->l;
- top += pow(i1 - i2, 2.0);
- bot += pow(i1, 2.0);
+ i1 = scale * lookup_intensity(ref1, h, k, l);
+ i2 = lookup_intensity(ref2, h, k, l);
- } /* else reflection not measured so don't worry about it */
+ top += pow(i1 - i2, 2.0);
+ bot += pow(i1, 2.0);
}
@@ -120,30 +134,29 @@ static double internal_r2(const double *ref1, const unsigned int *c1,
}
-static double internal_rmerge(const double *ref1, const unsigned int *c1,
- const double *ref2, const unsigned int *c2,
- double scale)
+static double internal_rmerge(const double *ref1, const double *ref2,
+ ReflItemList *items, double scale)
{
double top = 0.0;
double bot = 0.0;
int i;
- /* Start from i=1 -> skip central beam */
- for ( i=1; i<LIST_SIZE; i++ ) {
-
- if ( c1[i] && c2[i] ) {
+ for ( i=0; i<num_items(items); i++ ) {
- double f1, f2;
+ double i1, i2, f1, f2;
+ struct refl_item *it;
+ signed int h, k, l;
- if ( (ref1[i]<0.0) || (ref2[i]<0.0) ) continue;
+ it = get_item(items, i);
+ h = it->h; k = it->k; l = it->l;
- f1 = sqrt(ref1[i]) / (scale*(double)c1[i]);
- f2 = sqrt(ref2[i]) / (double)c2[i];
-
- top += fabs(f1 - f2);
- bot += f1 + f2;
+ 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;
- } /* else reflection not measured so don't worry about it */
+ top += fabs(f1 - f2);
+ bot += f1 + f2;
}
@@ -157,11 +170,9 @@ static double calc_r(double scale, void *params)
switch ( rp->fom) {
case R_MERGE :
- return internal_rmerge(rp->ref1, rp->c1,
- rp->ref2, rp->c2, scale);
+ return internal_rmerge(rp->ref1, rp->ref2, rp->items, scale);
case R_2 :
- return internal_r2(rp->ref1, rp->c1,
- rp->ref2, rp->c2, scale);
+ return internal_r2(rp->ref1, rp->ref2, rp->items, scale);
}
ERROR("No such FoM!\n");
@@ -169,9 +180,8 @@ static double calc_r(double scale, void *params)
}
-static double r_minimised(const double *ref1, const unsigned int *c1,
- const double *ref2, const unsigned int *c2,
- double *scalep, int fom)
+static double r_minimised(const double *ref1, const double *ref2,
+ ReflItemList *items, double *scalep, int fom)
{
gsl_function F;
gsl_min_fminimizer *s;
@@ -182,8 +192,7 @@ static double r_minimised(const double *ref1, const unsigned int *c1,
rp.ref1 = ref1;
rp.ref2 = ref2;
- rp.c1 = c1;
- rp.c2 = c2;
+ rp.items = items;
rp.fom = fom;
F.function = &calc_r;
@@ -194,10 +203,10 @@ static double r_minimised(const double *ref1, const unsigned int *c1,
/* Initial guess */
switch ( fom ) {
case R_MERGE :
- scale = stat_scale_sqrti(ref1, c1, ref2, c2);
+ scale = stat_scale_sqrti(ref1, ref2, items);
break;
case R_2 :
- scale = stat_scale_intensity(ref1, c1, ref2, c2);
+ scale = stat_scale_intensity(ref1, ref2, items);
break;
}
//STATUS("Initial scale factor estimate: %5.2e\n", scale);
@@ -236,52 +245,52 @@ static double r_minimised(const double *ref1, const unsigned int *c1,
}
-double stat_rmerge(const double *ref1, const unsigned int *c1,
- const double *ref2, const unsigned int *c2,
- double *scalep)
+double stat_rmerge(const double *ref1, const double *ref2,
+ ReflItemList *items, double *scalep)
{
- return r_minimised(ref1, c1, ref2, c2, scalep, R_MERGE);
+ return r_minimised(ref1, ref2, items, scalep, R_MERGE);
}
-double stat_r2(const double *ref1, const unsigned int *c1,
- const double *ref2, const unsigned int *c2,
- double *scalep)
+double stat_r2(const double *ref1, const double *ref2,
+ ReflItemList *items, double *scalep)
{
- return r_minimised(ref1, c1, ref2, c2, scalep, R_2);
+ return r_minimised(ref1, ref2, items, scalep, R_2);
}
-double stat_pearson(const double *ref1, const unsigned int *c1,
- const double *ref2, const unsigned int *c2)
+double stat_pearson(const double *ref1, const double *ref2, ReflItemList *items)
{
- double vec1[4096];
- double vec2[4096];
- signed int h, k, l;
+ double *vec1, *vec2;
int i = 0;
+ int ni = num_items(items);
+ double val;
- for ( l=-INDMAX; l<INDMAX; l++ ) {
- for ( k=-INDMAX; k<INDMAX; k++ ) {
- for ( h=-INDMAX; h<INDMAX; h++ ) {
+ vec1 = malloc(ni*sizeof(double));
+ vec2 = malloc(ni*sizeof(double));
- double i1, i2;
- unsigned int c1s, c2s;
+ for ( i=0; i<ni; i++ ) {
+
+ double i1, i2, f1, f2;
+ struct refl_item *it;
+ signed int h, k, l;
- c1s = lookup_count(c1, h, k, l);
- c2s = lookup_count(c2, h, k, l);
+ it = get_item(items, i);
+ h = it->h; k = it->k; l = it->l;
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;
- if ( c1s && c2s && (i1>0.0) && (i2>0.0) ) {
- vec1[i] = sqrt(i1 / (double)c1s);
- vec2[i] = sqrt(i2 / (double)c2s);
- i++;
- }
+ vec1[i] = f1;
+ vec2[i] = f2;
}
- }
- }
- return gsl_stats_correlation(vec1, 1, vec2, 1, i);
+ val = gsl_stats_correlation(vec1, 1, vec2, 1, i);
+ free(vec1);
+ free(vec2);
+
+ return val;
}
diff --git a/src/statistics.h b/src/statistics.h
index 5e3ca52c..c6c2f09c 100644
--- a/src/statistics.h
+++ b/src/statistics.h
@@ -16,18 +16,20 @@
#ifndef STATISTICS_H
#define STATISTICS_H
-extern double stat_scale_intensity(const double *ref1, const unsigned int *c1,
- const double *ref2, const unsigned int *c2);
-extern double stat_r2(const double *ref1, const unsigned int *c1,
- const double *ref2, const unsigned int *c2,
- double *scalep);
+#include "utils.h"
-extern double stat_rmerge(const double *ref1, const unsigned int *c1,
- const double *ref2, const unsigned int *c2,
- double *scalep);
+extern double stat_scale_intensity(const double *ref1, const double *ref2,
+ ReflItemList *items);
+
+extern double stat_rmerge(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_pearson(const double *ref1, const double *ref2,
+ ReflItemList *items);
-extern double stat_pearson(const double *ref1, const unsigned int *c1,
- const double *ref2, const unsigned int *c2);
#endif /* STATISTICS_H */
diff --git a/src/utils.c b/src/utils.c
index 21ba882f..0f20435e 100644
--- a/src/utils.c
+++ b/src/utils.c
@@ -387,24 +387,23 @@ int num_items(const ReflItemList *items)
}
-unsigned int *items_to_counts(ReflItemList *items)
+void union_op_items(ReflItemList *items, ReflItemList *newi)
{
- unsigned int *c;
- int i;
+ int n, i;
- c = new_list_count();
+ n = num_items(newi);
+ for ( i=0; i<n; i++ ) {
- for ( i=0; i<num_items(items); i++ ) {
- struct refl_item *r;
- r = get_item(items, i);
- set_count(c, r->h, r->k, r->l, 1);
- }
+ struct refl_item *r = get_item(newi, i);
+ if ( find_op(items, r->op) ) continue;
+
+ add_item_with_op(items, r->h, r->k, r->l, r->op);
- return c;
+ }
}
-void union_op_items(ReflItemList *items, ReflItemList *newi)
+void union_items(ReflItemList *items, ReflItemList *newi)
{
int n, i;
@@ -412,9 +411,28 @@ void union_op_items(ReflItemList *items, ReflItemList *newi)
for ( i=0; i<n; i++ ) {
struct refl_item *r = get_item(newi, i);
- if ( find_op(items, r->op) ) continue;
+ if ( find_item(items, r->h, r->k, r->l) ) continue;
add_item_with_op(items, r->h, r->k, r->l, r->op);
}
}
+
+
+ReflItemList *intersection_items(ReflItemList *i1, ReflItemList *i2)
+{
+ int n, i;
+ ReflItemList *res = new_items();
+
+ n = num_items(i1);
+ for ( i=0; i<n; i++ ) {
+
+ struct refl_item *r = get_item(i1, i);
+ if ( find_item(i2, r->h, r->k, r->l) ) {
+ add_item_with_op(res, r->h, r->k, r->l, r->op);
+ }
+
+ }
+
+ return res;
+}
diff --git a/src/utils.h b/src/utils.h
index a7433d28..f7e98469 100644
--- a/src/utils.h
+++ b/src/utils.h
@@ -200,6 +200,8 @@ extern struct refl_item *get_item(ReflItemList *items, int i);
extern int num_items(const ReflItemList *items);
extern unsigned int *items_to_counts(ReflItemList *items);
extern void union_op_items(ReflItemList *items, ReflItemList *newi);
+extern void union_items(ReflItemList *items, ReflItemList *newi);
+extern ReflItemList *intersection_items(ReflItemList *i1, ReflItemList *i2);
/* ------------------------------ Message macros ---------------------------- */
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-15 22:27:53 +0000