aboutsummaryrefslogtreecommitdiff
path: root/src/rejection.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2018-08-20 16:42:25 +0200
committerThomas White <taw@physics.org>2018-08-30 17:18:53 +0200
commit609dcac27949236c7a1acda0bd3e4342cdb328f6 (patch)
treea87bcfc5fbfcbe4dcdd578ef95cc7c9d35cbee53 /src/rejection.c
parent92f2ebfd11d0ccd91ab299f95ffb2d99a457643d (diff)
Initial delta CC half stuff
Diffstat (limited to 'src/rejection.c')
-rw-r--r--src/rejection.c211
1 files changed, 165 insertions, 46 deletions
diff --git a/src/rejection.c b/src/rejection.c
index 9d41f9b4..716baa1b 100644
--- a/src/rejection.c
+++ b/src/rejection.c
@@ -94,74 +94,194 @@ void early_rejection(Crystal **crystals, int n)
}
-static void check_cc(Crystal *cr, RefList *full)
+struct contribs
+{
+ int serial; /* h,k,l */
+ int n_contrib;
+ int max_contrib;
+ Reflection **contribs;
+ Crystal **contrib_crystals;
+};
+
+
+struct contributionlist
+{
+ struct contribs *contribs;
+ int n; /* Number of reflections */
+};
+
+
+static int alloc_contribs(struct contribs *c)
+{
+ c->contribs = realloc(c->contribs, c->max_contrib*sizeof(Reflection *));
+ c->contrib_crystals = realloc(c->contrib_crystals,
+ c->max_contrib*sizeof(Crystal *));
+ if ( c->contribs == NULL ) return 1;
+ if ( c->contrib_crystals == NULL ) return 1;
+ return 0;
+}
+
+
+static void free_contributions(struct contributionlist *clist)
+{
+ int i;
+ for ( i=0; i<clist->n; i++ ) {
+ free(clist->contribs[i].contribs);
+ free(clist->contribs[i].contrib_crystals);
+ }
+ free(clist->contribs);
+ free(clist);
+}
+
+
+static struct contributionlist *find_all_contributions(Crystal **crystals,
+ int n, RefList *full)
{
- RefList *list = crystal_get_reflections(cr);
Reflection *refl;
RefListIterator *iter;
- double G = crystal_get_osf(cr);
- double B = crystal_get_Bfac(cr);
- UnitCell *cell = crystal_get_cell(cr);
- double *vec1, *vec2;
- int n, i;
- double cc;
-
- n = num_reflections(list);
- vec1 = malloc(n*sizeof(double));
- vec2 = malloc(n*sizeof(double));
- if ( (vec1 == NULL) || (vec2 == NULL) ) {
- ERROR("Not enough memory to check CCs\n");
- return;
+ struct contributionlist *clist;
+ int nc = 0;
+
+ clist = malloc(sizeof(struct contributionlist));
+ if ( clist == NULL ) return NULL;
+
+ clist->n = num_reflections(full);
+ clist->contribs = malloc(clist->n*sizeof(struct contribs));
+ if ( clist->contribs == NULL ) {
+ free(clist);
+ return NULL;
}
- i = 0;
- for ( refl = first_refl(list, &iter);
+ for ( refl = first_refl(full, &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
+ int i;
signed int h, k, l;
- double pobs, pcalc;
- double res, corr, Ipart;
- Reflection *match;
-
- if ( !get_flag(refl) ) continue; /* Not free-flagged */
+ struct contribs *c;
get_indices(refl, &h, &k, &l);
- res = resolution(cell, h, k, l);
- if ( 2.0*res > crystal_get_resolution_limit(cr) ) continue;
+ c = &clist->contribs[nc++];
+
+ c->serial = SERIAL(h, k, l);
+ c->n_contrib = 0;
+ c->max_contrib = 32;
+ c->contribs = NULL;
+ c->contrib_crystals = NULL;
+ if ( alloc_contribs(c) ) return NULL;
- match = find_refl(full, h, k, l);
- if ( match == NULL ) continue;
+ /* Find all versions of this reflection */
+ for ( i=0; i<n; i++ ) {
- /* Calculated partiality */
- pcalc = get_partiality(refl);
+ Reflection *refl2;
+ RefList *list2 = crystal_get_reflections(crystals[i]);
+ refl2 = find_refl(list2, h, k, l);
- /* Observed partiality */
- corr = G * exp(B*res*res) * get_lorentz(refl);
- Ipart = get_intensity(refl) * corr;
- pobs = Ipart / get_intensity(match);
+ if ( refl2 == NULL ) continue;
+ do {
+ c->contribs[c->n_contrib] = refl2;
+ c->contrib_crystals[c->n_contrib] = crystals[i];
+ c->n_contrib++;
- vec1[i] = pobs;
- vec2[i] = pcalc;
- i++;
+ if ( c->n_contrib == c->max_contrib ) {
+ c->max_contrib += 64;
+ if ( alloc_contribs(c) ) return NULL;
+ }
+
+ refl2 = next_found_refl(refl2);
+
+ } while ( refl2 != NULL );
+ }
+
+ progress_bar(nc, clist->n, "Collecting contributions");
}
- cc = gsl_stats_correlation(vec1, 1, vec2, 1, i);
- //printf("%f\n", cc);
- if ( cc < 0.5 ) crystal_set_user_flag(cr, PRFLAG_CC);
+ return clist;
+}
+
- free(vec1);
- free(vec2);
+static struct contribs *lookup_contribs(struct contributionlist *clist,
+ signed int h, signed int k, signed int l)
+{
+ int i, serial;
+
+ serial = SERIAL(h, k, l);
+
+ for ( i=0; i<clist->n; i++ ) {
+ if ( clist->contribs[i].serial == serial ) {
+ return &clist->contribs[i];
+ }
+ }
+ return NULL;
}
-static UNUSED void check_ccs(Crystal **crystals, int n_crystals, RefList *full)
+static double calculate_cchalf(Crystal **crystals, int n,
+ struct contributionlist *clist, Crystal *exclude)
{
int i;
+ double acc_var = 0.0;
+ double acc_mean = 0.0;
+ int acc_n = 0;
+ double sigE, sigy;
+
+ for ( i=0; i<n; i++ ) {
+
+ Crystal *cr = crystals[i];
+ RefList *list = crystal_get_reflections(cr);
+ Reflection *refl;
+ RefListIterator *iter;
+
+ if ( cr == exclude ) continue;
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ signed int h, k, l;
+ struct contribs *c;
+ int j;
+
+ get_indices(refl, &h, &k, &l);
+ c = lookup_contribs(clist, h, k, l);
+ assert(c != NULL);
+ for ( j=0; j<c->n_contrib; j++ ) {
+
+ Reflection *reflc;
+ double Ii;
+
+ if ( c->contrib_crystals[j] == exclude ) {
+ continue;
+ }
+
+ reflc = c->contribs[j];
+
+ /* FIXME: Apply corrections */
+ Ii = get_intensity(reflc);
+ acc_mean += Ii;
+ acc_var += Ii * Ii;
+ acc_n++;
+ }
+
+ }
- for ( i=0; i<n_crystals; i++ ) {
- check_cc(crystals[i], full);
}
+
+
+ return 0.0; /* FIXME */
+}
+
+
+static void check_deltacchalf(Crystal **crystals, int n, RefList *full)
+{
+ struct contributionlist *clist;
+
+ clist = find_all_contributions(crystals, n, full);
+
+ STATUS("Overall CChalf = %f\n",
+ calculate_cchalf(crystals, n, clist, NULL));
+
+ free_contributions(clist);
}
@@ -197,9 +317,8 @@ void check_rejection(Crystal **crystals, int n, RefList *full, double max_B)
int i;
int n_acc = 0;
- /* Check according to CCs FIXME: Disabled */
- //if ( full != NULL ) check_ccs(crystals, n, full);
-
+ /* Check according to delta CC½ */
+ if ( full != NULL ) check_deltacchalf(crystals, n, full);
for ( i=0; i<n; i++ ) {
if ( fabs(crystal_get_Bfac(crystals[i])) > max_B ) {