aboutsummaryrefslogtreecommitdiff
path: root/src/get_hkl.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/get_hkl.c')
-rw-r--r--src/get_hkl.c283
1 files changed, 145 insertions, 138 deletions
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;
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-14 01:16:25 +0000