aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/merge.c55
-rw-r--r--src/merge.h6
-rw-r--r--src/partialator.c60
-rw-r--r--src/post-refinement.c5
-rw-r--r--src/rejection.c27
5 files changed, 102 insertions, 51 deletions
diff --git a/src/merge.c b/src/merge.c
index 56ca2020..8017182b 100644
--- a/src/merge.c
+++ b/src/merge.c
@@ -187,7 +187,7 @@ static void run_merge_job(void *vwargs, int cookie)
Reflection *f;
signed int h, k, l;
double mean, sumweight, M2, temp, delta, R;
- double corr, res, w;
+ double res, w;
struct reflection_contributions *c;
if ( get_partiality(refl) < MIN_PART_MERGE ) continue;
@@ -217,25 +217,16 @@ static void run_merge_job(void *vwargs, int cookie)
continue;
}
- /* Total (multiplicative) correction factor */
- corr = 1.0/G * exp(B*res*res) * get_lorentz(refl)
- / get_partiality(refl);
- if ( isnan(corr) ) {
- ERROR("NaN corr:\n");
- ERROR("G = %f, B = %e\n", G, B);
- ERROR("res = %e\n", res);
- ERROR("p = %f\n", get_partiality(refl));
- unlock_reflection(f);
- continue;
- }
-
/* Reflections count less the more they have to be scaled up */
- w = 1.0;
+ w = get_partiality(refl);
/* Running mean and variance calculation */
temp = w + sumweight;
- if (ln_merge) delta = log(get_intensity(refl)*corr) - mean;
- else delta = get_intensity(refl)*corr - mean;
+ if ( ln_merge ) {
+ delta = log(correct_reflection(refl, G, B, res)) - mean;
+ } else {
+ delta = correct_reflection(refl, G, B, res) - mean;
+ }
R = delta * w / temp;
set_intensity(f, mean + R);
set_temp2(f, M2 + sumweight * delta * R);
@@ -352,6 +343,25 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
}
+/* Correct intensity in pattern for scaling and Lorentz factors,
+ * but not partiality nor polarisation */
+double correct_reflection_nopart(Reflection *refl, double osf, double Bfac,
+ double res)
+{
+ double corr = osf * exp(-Bfac*res*res);
+ return (get_intensity(refl) / corr) / get_lorentz(refl);
+}
+
+
+/* Correct intensity in pattern for scaling, partiality and Lorentz factors,
+ * but not polarisation */
+double correct_reflection(Reflection *refl, double osf, double Bfac, double res)
+{
+ double Ipart = correct_reflection_nopart(refl, osf, Bfac, res);
+ return Ipart / get_partiality(refl);
+}
+
+
double residual(Crystal *cr, const RefList *full, int free,
int *pn_used, const char *filename)
{
@@ -368,7 +378,7 @@ double residual(Crystal *cr, const RefList *full, int free,
refl != NULL;
refl = next_refl(refl, iter) )
{
- double p, w, corr, res;
+ double w, res;
signed int h, k, l;
Reflection *match;
double I_full;
@@ -384,13 +394,9 @@ double residual(Crystal *cr, const RefList *full, int free,
if ( get_redundancy(match) < 2 ) continue;
- p = get_partiality(refl);
- //if ( p < 0.2 ) continue;
-
- corr = get_lorentz(refl) / (G * exp(-B*res*res));
- int1 = get_intensity(refl) * corr;
- int2 = p*I_full;
- w = p;
+ int1 = correct_reflection_nopart(refl, G, B, res);
+ int2 = get_partiality(refl)*I_full;
+ w = 1.0;
num += fabs(int1 - int2) * w;
den += fabs(int1 + int2) * w/2.0;
@@ -454,6 +460,7 @@ double log_residual(Crystal *cr, const RefList *full, int free,
fx = log(G) - B*s*s + log(p) + log(I_full) - log(I_partial) - log(L);
w = 1.0;
dev += w*fx*fx;
+ n_used++;
if ( fh != NULL ) {
fprintf(fh, "%4i %4i %4i %e %e\n",
diff --git a/src/merge.h b/src/merge.h
index 476f8755..2c2cfeda 100644
--- a/src/merge.h
+++ b/src/merge.h
@@ -47,6 +47,12 @@ extern RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
int min_meas, double push_res, int use_weak,
int ln_merge);
+extern double correct_reflection_nopart(Reflection *refl, double osf,
+ double Bfac, double res);
+
+extern double correct_reflection(Reflection *refl, double osf, double Bfac,
+ double res);
+
extern double residual(Crystal *cr, const RefList *full, int free,
int *pn_used, const char *filename);
diff --git a/src/partialator.c b/src/partialator.c
index 86431e4c..8cf1dad9 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -638,7 +638,7 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
{
signed int h, k, l;
double pobs, pcalc;
- double res, corr, Ipart;
+ double res, Ipart;
Reflection *match;
if ( !get_flag(refl) ) continue; /* Not free-flagged */
@@ -660,8 +660,7 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
pcalc = get_partiality(refl);
/* Observed partiality */
- corr = G * exp(B*res*res) * get_lorentz(refl);
- Ipart = get_intensity(refl) * corr;
+ Ipart = correct_reflection_nopart(refl, G, B, res);
pobs = Ipart / get_intensity(match);
fprintf(fh, "%5i %4i %4i %4i %e %e %8.3f %8.3f %s\n",
@@ -716,6 +715,10 @@ static void all_residuals(Crystal **crystals, int n_crystals, RefList *full,
int n_nan_linear_free = 0;
int n_nan_log = 0;
int n_nan_log_free = 0;
+ int n_non_linear = 0;
+ int n_non_linear_free = 0;
+ int n_non_log = 0;
+ int n_non_log_free = 0;
*presidual = 0.0;
*pfree_residual = 0.0;
@@ -725,19 +728,35 @@ static void all_residuals(Crystal **crystals, int n_crystals, RefList *full,
for ( i=0; i<n_crystals; i++ ) {
double r, free_r, log_r, free_log_r;
+ int n;
if ( crystal_get_user_flag(crystals[i]) ) continue;
/* Scaling should have been done right before calling this */
- r = residual(crystals[i], full, 0, NULL, NULL);
- free_r = residual(crystals[i], full, 1, NULL, NULL);
- log_r = log_residual(crystals[i], full, 0, NULL, NULL);
- free_log_r = log_residual(crystals[i], full, 1, NULL, NULL);
-
- if ( isnan(r) ) n_nan_linear++;
- if ( isnan(free_r) ) n_nan_linear_free++;
- if ( isnan(log_r) ) n_nan_log++;
- if ( isnan(free_log_r) ) n_nan_log_free++;
+ r = residual(crystals[i], full, 0, &n, NULL);
+ if ( n == 0 ) {
+ n_non_linear++;
+ } else if ( isnan(r) ) {
+ n_nan_linear++;
+ }
+ free_r = residual(crystals[i], full, 1, &n, NULL);
+ if ( n == 0 ) {
+ n_non_linear_free++;
+ } else if ( isnan(free_r) ) {
+ n_nan_linear_free++;
+ }
+ log_r = log_residual(crystals[i], full, 0, &n, NULL);
+ if ( n == 0 ) {
+ n_non_log++;
+ } else if ( isnan(log_r) ) {
+ n_nan_log++;
+ }
+ free_log_r = log_residual(crystals[i], full, 1, &n, NULL);
+ if ( n == 0 ) {
+ n_non_log_free++;
+ } else if ( isnan(free_log_r) ) {
+ n_nan_log_free++;
+ }
if ( isnan(r) || isnan(free_r)
|| isnan(log_r) || isnan(free_log_r) ) continue;
@@ -750,6 +769,23 @@ static void all_residuals(Crystal **crystals, int n_crystals, RefList *full,
n_used++;
}
+ if ( n_non_linear ) {
+ ERROR("WARNING: %i crystals had no reflections in linear "
+ "residual calculation\n", n_non_linear);
+ }
+ if ( n_non_linear_free ) {
+ ERROR("WARNING: %i crystals had no reflections in linear free "
+ "residual calculation\n", n_non_linear_free);
+ }
+ if ( n_non_log ) {
+ ERROR("WARNING: %i crystals had no reflections in log "
+ "residual calculation\n", n_non_log);
+ }
+ if ( n_non_log_free ) {
+ ERROR("WARNING: %i crystals had no reflections in log free "
+ "residual calculation\n", n_non_log_free);
+ }
+
if ( n_nan_linear ) {
ERROR("WARNING: %i crystals had NaN linear residuals\n",
n_nan_linear);
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 0a799b9e..b1af92db 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -559,7 +559,7 @@ void write_specgraph(Crystal *crystal, const RefList *full,
refl != NULL;
refl = next_refl(refl, iter) )
{
- double corr, Ipart, Ifull, pobs, pcalc;
+ double Ipart, Ifull, pobs, pcalc;
double res;
signed int h, k, l;
Reflection *match;
@@ -573,8 +573,7 @@ void write_specgraph(Crystal *crystal, const RefList *full,
/* Don't calculate pobs if reference reflection is weak */
if ( fabs(get_intensity(match)) / get_esd_intensity(match) < 3.0 ) continue;
- corr = G * exp(B*res*res) * get_lorentz(refl);
- Ipart = get_intensity(refl) * corr;
+ Ipart = correct_reflection_nopart(refl, G, B, res);
Ifull = get_intensity(match);
pobs = Ipart / Ifull;
pcalc = get_partiality(refl);
diff --git a/src/rejection.c b/src/rejection.c
index b3a3cb67..5d1041fa 100644
--- a/src/rejection.c
+++ b/src/rejection.c
@@ -143,11 +143,7 @@ static int calculate_refl_mean_var(RefList *full)
G = crystal_get_osf(c->contrib_crystals[j]);
B = crystal_get_Bfac(c->contrib_crystals[j]);
- Ii = get_intensity(c->contribs[j]);
-
- /* Total (multiplicative) correction factor */
- Ii *= 1.0/G * exp(B*res*res) * get_lorentz(c->contribs[j])
- / get_partiality(c->contribs[j]);
+ Ii = correct_reflection(c->contribs[j], G, B, res);
Ex += Ii - K;
Ex2 += (Ii - K) * (Ii - K);
@@ -232,16 +228,13 @@ static double calculate_cchalf(RefList *template, RefList *full,
while ( exrefl != NULL ) {
double G, B;
- double Ii = get_intensity(exrefl);
G = crystal_get_osf(exclude);
B = crystal_get_Bfac(exclude);
if ( get_partiality(exrefl) > MIN_PART_MERGE ) {
- /* Total (multiplicative) correction factor */
- Ii *= 1.0/G * exp(B*res*res) * get_lorentz(exrefl)
- / get_partiality(exrefl);
+ double Ii = correct_reflection(exrefl, G, B, res);
/* Remove contribution of this reflection */
Ex -= Ii - K;
@@ -289,6 +282,7 @@ static void check_deltacchalf(Crystal **crystals, int n, RefList *full)
double mean, sd;
int nref = 0;
int nnan = 0;
+ int nnon = 0;
if ( calculate_refl_mean_var(full) ) {
STATUS("No reflection contributions for deltaCChalf "
@@ -315,13 +309,22 @@ static void check_deltacchalf(Crystal **crystals, int n, RefList *full)
//STATUS("Without = %f", cchalfi*100.0);
//STATUS(" Delta = %f ", (cchalf - cchalfi)*100.0);
//STATUS("(nref = %i)\n", nref);
- vals[i] = cchalf - cchalfi;
- if ( isnan(vals[i]) || isinf(vals[i]) ) {
+ if ( nref == 0 ) {
vals[i] = 0.0;
- nnan++;
+ nnon++;
+ } else {
+ vals[i] = cchalf - cchalfi;
+ if ( isnan(vals[i]) || isinf(vals[i]) ) {
+ vals[i] = 0.0;
+ nnan++;
+ }
}
progress_bar(i, n-1, "Calculating deltaCChalf");
}
+ if ( nnon > 0 ) {
+ STATUS("WARNING: %i patterns had no reflections in deltaCChalf "
+ "calculation (I set deltaCChalf=zero for them)\n", nnon);
+ }
if ( nnan > 0 ) {
STATUS("WARNING: %i NaN or inf deltaCChalf values were "
"replaced with zero\n", nnan);
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-04 01:52:57 +0000