From 5790b06b2e0080c48e1e9a33eb0b43914f2b5824 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 2 May 2018 17:41:47 +0200 Subject: Move residual() and log_residual() to merge.c --- src/merge.c | 118 ++++++++++++++++++++++++++++++++++++++++++++++++++ src/merge.h | 6 +++ src/post-refinement.c | 53 +---------------------- src/post-refinement.h | 3 -- src/scaling.c | 66 ---------------------------- src/scaling.h | 3 -- 6 files changed, 125 insertions(+), 124 deletions(-) diff --git a/src/merge.c b/src/merge.c index 863492f6..aab47ee6 100644 --- a/src/merge.c +++ b/src/merge.c @@ -282,3 +282,121 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads, reflist_free(full); return full2; } + + +double residual(Crystal *cr, const RefList *full, int free, + int *pn_used, const char *filename) +{ + Reflection *refl; + RefListIterator *iter; + int n_used = 0; + double num = 0.0; + double den = 0.0; + double G = crystal_get_osf(cr); + double B = crystal_get_Bfac(cr); + UnitCell *cell = crystal_get_cell(cr); + + for ( refl = first_refl(crystal_get_reflections(cr), &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + double p, w, corr, res; + signed int h, k, l; + Reflection *match; + double I_full; + double int1, int2; + + if ( free && !get_flag(refl) ) continue; + + get_indices(refl, &h, &k, &l); + res = resolution(cell, h, k, l); + match = find_refl(full, h, k, l); + if ( match == NULL ) continue; + I_full = get_intensity(match); + + if ( get_redundancy(match) < 2 ) continue; + + p = get_partiality(refl); + //if ( p < 0.2 ) continue; + + corr = G * exp(B*res*res) * get_lorentz(refl); + int1 = get_intensity(refl) * corr; + int2 = p*I_full; + w = p; + + num += fabs(int1 - int2) * w; + den += fabs(int1 + int2) * w/2.0; + + n_used++; + + } + + if ( pn_used != NULL ) *pn_used = n_used; + return num/(den*sqrt(2)); +} + + +double log_residual(Crystal *cr, const RefList *full, int free, + int *pn_used, const char *filename) +{ + double dev = 0.0; + double G, B; + Reflection *refl; + RefListIterator *iter; + int n_used = 0; + FILE *fh = NULL; + + G = crystal_get_osf(cr); + B = crystal_get_Bfac(cr); + if ( filename != NULL ) { + fh = fopen(filename, "a"); + if ( fh == NULL ) { + ERROR("Failed to open '%s'\n", filename); + } + } + + for ( refl = first_refl(crystal_get_reflections(cr), &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + double p, L, s, w; + signed int h, k, l; + Reflection *match; + double esd, I_full, I_partial; + double fx, dc; + + if ( free && !get_flag(refl) ) continue; + + get_indices(refl, &h, &k, &l); + match = find_refl(full, h, k, l); + if ( match == NULL ) continue; + + p = get_partiality(refl); + L = get_lorentz(refl); + I_partial = get_intensity(refl); + I_full = get_intensity(match); + esd = get_esd_intensity(refl); + s = resolution(crystal_get_cell(cr), h, k, l); + + if ( I_partial <= 3.0*esd ) continue; /* Also because of log */ + if ( get_redundancy(match) < 2 ) continue; + if ( I_full <= 0 ) continue; /* Because log */ + if ( p <= 0.0 ) continue; /* Because of log */ + + fx = -log(G) + log(p) - log(L) - B*s*s + log(I_full); + dc = log(I_partial) - fx; + w = 1.0; + dev += w*dc*dc; + + if ( fh != NULL ) { + fprintf(fh, "%4i %4i %4i %e %e\n", + h, k, l, s, dev); + } + + } + + if ( fh != NULL ) fclose(fh); + + if ( pn_used != NULL ) *pn_used = n_used; + return dev; +} diff --git a/src/merge.h b/src/merge.h index 56229e4c..b33afdea 100644 --- a/src/merge.h +++ b/src/merge.h @@ -42,4 +42,10 @@ extern RefList *merge_intensities(Crystal **crystals, int n, int n_threads, int min_meas, double push_res, int use_weak); +extern double residual(Crystal *cr, const RefList *full, int free, + int *pn_used, const char *filename); + +extern double log_residual(Crystal *cr, const RefList *full, int free, + int *pn_used, const char *filename); + #endif /* MERGE */ diff --git a/src/post-refinement.c b/src/post-refinement.c index 293fe332..4e83eacc 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -45,6 +45,7 @@ #include "cell-utils.h" #include "reflist-utils.h" #include "scaling.h" +#include "merge.h" struct prdata @@ -83,58 +84,6 @@ const char *str_prflag(enum prflag flag) } -double residual(Crystal *cr, const RefList *full, int free, - int *pn_used, const char *filename) -{ - Reflection *refl; - RefListIterator *iter; - int n_used = 0; - double num = 0.0; - double den = 0.0; - double G = crystal_get_osf(cr); - double B = crystal_get_Bfac(cr); - UnitCell *cell = crystal_get_cell(cr); - - for ( refl = first_refl(crystal_get_reflections(cr), &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - double p, w, corr, res; - signed int h, k, l; - Reflection *match; - double I_full; - double int1, int2; - - if ( free && !get_flag(refl) ) continue; - - get_indices(refl, &h, &k, &l); - res = resolution(cell, h, k, l); - match = find_refl(full, h, k, l); - if ( match == NULL ) continue; - I_full = get_intensity(match); - - if ( get_redundancy(match) < 2 ) continue; - - p = get_partiality(refl); - //if ( p < 0.2 ) continue; - - corr = G * exp(B*res*res) * get_lorentz(refl); - int1 = get_intensity(refl) * corr; - int2 = p*I_full; - w = p; - - num += fabs(int1 - int2) * w; - den += fabs(int1 + int2) * w/2.0; - - n_used++; - - } - - if ( pn_used != NULL ) *pn_used = n_used; - return num/(den*sqrt(2)); -} - - static UnitCell *rotate_cell_xy(const UnitCell *cell, double ang1, double ang2) { UnitCell *o; diff --git a/src/post-refinement.h b/src/post-refinement.h index b8923d2c..71a6d7f3 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -73,9 +73,6 @@ extern void write_specgraph(Crystal *crystal, const RefList *full, extern double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel); -extern double residual(Crystal *cr, const RefList *full, int free, - int *pn_used, const char *filename); - extern void write_test_logs(Crystal *crystal, const RefList *full, signed int cycle, int serial); diff --git a/src/scaling.c b/src/scaling.c index 5f896e65..92f618ae 100644 --- a/src/scaling.c +++ b/src/scaling.c @@ -221,72 +221,6 @@ static double scale_iterate(Crystal *cr, const RefList *full, int *nr) } -double log_residual(Crystal *cr, const RefList *full, int free, - int *pn_used, const char *filename) -{ - double dev = 0.0; - double G, B; - Reflection *refl; - RefListIterator *iter; - int n_used = 0; - FILE *fh = NULL; - - G = crystal_get_osf(cr); - B = crystal_get_Bfac(cr); - if ( filename != NULL ) { - fh = fopen(filename, "a"); - if ( fh == NULL ) { - ERROR("Failed to open '%s'\n", filename); - } - } - - for ( refl = first_refl(crystal_get_reflections(cr), &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - double p, L, s, w; - signed int h, k, l; - Reflection *match; - double esd, I_full, I_partial; - double fx, dc; - - if ( free && !get_flag(refl) ) continue; - - get_indices(refl, &h, &k, &l); - match = find_refl(full, h, k, l); - if ( match == NULL ) continue; - - p = get_partiality(refl); - L = get_lorentz(refl); - I_partial = get_intensity(refl); - I_full = get_intensity(match); - esd = get_esd_intensity(refl); - s = resolution(crystal_get_cell(cr), h, k, l); - - if ( I_partial <= 3.0*esd ) continue; /* Also because of log */ - if ( get_redundancy(match) < 2 ) continue; - if ( I_full <= 0 ) continue; /* Because log */ - if ( p <= 0.0 ) continue; /* Because of log */ - - fx = -log(G) + log(p) - log(L) - B*s*s + log(I_full); - dc = log(I_partial) - fx; - w = 1.0; - dev += w*dc*dc; - - if ( fh != NULL ) { - fprintf(fh, "%4i %4i %4i %e %e\n", - h, k, l, s, dev); - } - - } - - if ( fh != NULL ) fclose(fh); - - if ( pn_used != NULL ) *pn_used = n_used; - return dev; -} - - static void do_scale_refine(Crystal *cr, const RefList *full, int *nr) { int i, done; diff --git a/src/scaling.h b/src/scaling.h index 25aca31e..b6958b6c 100644 --- a/src/scaling.h +++ b/src/scaling.h @@ -43,9 +43,6 @@ enum ScaleFlags SCALE_NO_B, /* Don't use Debye-Waller part */ }; -extern double log_residual(Crystal *cr, const RefList *full, int free, - int *pn_used, const char *filename); - extern int scale_one(const RefList *list1, const RefList *list2, int flags, double *G, double *B); -- cgit v1.2.3