aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2015-05-15 16:06:06 +0200
committerThomas White <taw@physics.org>2015-05-19 13:57:52 +0200
commit45475fd360e966ca85dc11fd8f15307073f767ab (patch)
treeb135c5c150e6f8f0967f3627c7aa02876d6dfdb0
parent6a9185eedebc6c2002b7d6ca2467d5e6a7f245f0 (diff)
More debugging / logarithm stuff
-rw-r--r--src/hrs-scaling.c8
-rw-r--r--src/partialator.c2
-rw-r--r--src/post-refinement.c89
3 files changed, 57 insertions, 42 deletions
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index 2a50d06a..f5067ea9 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -328,11 +328,11 @@ static void run_merge_job(void *vwargs, int cookie)
res = resolution(crystal_get_cell(cr), h, k, l);
/* Total (multiplicative) correction factor */
- corr = exp(2.0*B*res*res) * get_lorentz(refl)
- / (G * get_partiality(refl));
+ corr = exp(-G) * exp(B*res*res) * get_lorentz(refl)
+ / get_partiality(refl);
- esd = fabs(get_esd_intensity(refl) * corr);
- w = 1.0 / pow(esd, 0.5);
+ esd = get_esd_intensity(refl) * corr;
+ w = 1.0;
/* Running mean and variance calculation */
temp = w + sumweight;
diff --git a/src/partialator.c b/src/partialator.c
index b1e32c2f..7a6d3b05 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -273,7 +273,7 @@ static int set_initial_params(Crystal *cr, FILE *fh)
} else {
- crystal_set_osf(cr, 1.0);
+ crystal_set_osf(cr, 0.0);
crystal_set_Bfac(cr, 0.0);
}
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 3b584617..49eb12b1 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -192,7 +192,7 @@ double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
get_partial(refl, &rlow, &rhigh, &p);
if ( k == GPARAM_OSF ) {
- return 1.0/osf;
+ return 1.0;
}
if ( k == GPARAM_BFAC ) {
@@ -379,7 +379,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
refl = next_refl(refl, iter) )
{
signed int ha, ka, la;
- double I_full, delta_I;
+ double I_full, delta_I, esd;
double w;
double I_partial;
int k;
@@ -393,23 +393,22 @@ static double pr_iterate(Crystal *cr, const RefList *full,
match = find_refl(full, ha, ka, la);
if ( match == NULL ) continue;
- if ( (get_intensity(refl) < 3.0*get_esd_intensity(refl))
- || (get_partiality(refl) < MIN_PART_REFINE)
- || (get_redundancy(match) < 2) ) continue;
-
+ /* Merged intensitty */
I_full = get_intensity(match);
+ /* Actual measurement of this reflection from this pattern */
+ I_partial = get_intensity(refl);
+ esd = get_esd_intensity(refl);
+
+ if ( (get_partiality(refl) < MIN_PART_REFINE)
+ || (get_redundancy(match) < 2)
+ || (I_full <= 0) || (I_partial < 3*esd) ) continue;
+
p = get_partiality(refl);
L = get_lorentz(refl);
s = resolution(crystal_get_cell(cr), ha, ka, la);
- /* Actual measurement of this reflection from this pattern? */
- I_partial = get_intensity(refl);
-
/* Calculate the weight for this reflection */
- //w = pow(get_esd_intensity(refl), 2.0);
- //w += (p/L) * I_full * pow(get_esd_intensity(match), 2.0);
- //w = pow(w, -1.0);
w = 1.0;
/* Calculate all gradients for this reflection */
@@ -437,7 +436,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
}
- fx = log(G) + log(p) - log(L) - B*s*s + log(I_full);
+ fx = G + log(p) - log(L) - B*s*s + log(I_full);
delta_I = log(I_partial) - fx;
v_c = w * delta_I * gradients[k];
v_curr = gsl_vector_get(v, k);
@@ -483,49 +482,59 @@ static double pr_iterate(Crystal *cr, const RefList *full,
}
-static double guide_dev(Crystal *cr, const RefList *full)
+static double residual(Crystal *cr, const RefList *full, int verbose)
{
double dev = 0.0;
double G, B;
+ Reflection *refl;
+ RefListIterator *iter;
+ FILE *fh = NULL;
+
+ if ( verbose ) {
+ fh = fopen("residual.dat", "w");
+ }
G = crystal_get_osf(cr);
B = crystal_get_Bfac(cr);
- /* For each reflection */
- Reflection *refl;
- RefListIterator *iter;
-
for ( refl = first_refl(crystal_get_reflections(cr), &iter);
refl != NULL;
- refl = next_refl(refl, iter) ) {
-
- double p, L, s;
+ refl = next_refl(refl, iter) )
+ {
+ double p, L, s, w;
signed int h, k, l;
- Reflection *full_version;
- double I_full, I_partial;
+ Reflection *match;
+ double esd, I_full, I_partial;
double fx, dc;
- if ( (get_intensity(refl) < 3.0*get_esd_intensity(refl))
- || (get_partiality(refl) < MIN_PART_REFINE) ) continue;
-
get_indices(refl, &h, &k, &l);
- assert((h!=0) || (k!=0) || (l!=0));
-
- full_version = find_refl(full, h, k, l);
- if ( full_version == NULL ) continue;
+ 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(full_version);
+ I_full = get_intensity(match);
+ esd = get_esd_intensity(refl);
s = resolution(crystal_get_cell(cr), h, k, l);
- fx = log(G) + log(p) - log(L) - B*s*s + log(I_full);
+ if ( (get_partiality(refl) < MIN_PART_REFINE)
+ || (get_redundancy(match) < 2)
+ || (I_full <= 0) || (I_partial < 3*esd) ) continue;
+
+ fx = G + log(p) - log(L) - B*s*s + log(I_full);
dc = log(I_partial) - fx;
- dev += dc*dc;
+ w = 1.0;
+ dev += w*dc*dc;
+ if ( fh != NULL ) {
+ fprintf(fh, "%4i %4i %4i %e %.2f %e %f %f\n",
+ h, k, l, s, G, B, I_partial, I_full);
+ }
}
+ if ( fh != NULL ) fclose(fh);
+
return dev;
}
@@ -545,10 +554,11 @@ struct prdata pr_refine(Crystal *cr, const RefList *full,
if ( crystal_get_user_flag(cr) != 0 ) return prdata;
if ( verbose ) {
- dev = guide_dev(cr, full);
+ dev = residual(cr, full, 1);
STATUS("\n"); /* Deal with progress bar */
- STATUS("Before iteration: dev = %10.5e\n",
- dev);
+ STATUS("Initial G=%.2f, B=%e\n",
+ crystal_get_osf(cr), crystal_get_Bfac(cr));
+ STATUS("Initial dev = %10.5e\n", dev);
}
i = 0;
@@ -567,7 +577,7 @@ struct prdata pr_refine(Crystal *cr, const RefList *full,
update_partialities(cr, pmodel);
if ( verbose ) {
- dev = guide_dev(cr, full);
+ dev = residual(cr, full, 0);
STATUS("PR Iteration %2i: dev = %10.5e\n", i+1, dev);
}
@@ -579,5 +589,10 @@ struct prdata pr_refine(Crystal *cr, const RefList *full,
prdata.refined = 1;
}
+ if ( verbose ) {
+ STATUS("Final G=%.2f, B=%e\n", crystal_get_osf(cr),
+ crystal_get_Bfac(cr));
+ }
+
return prdata;
}