diff options
author | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-09-07 22:32:35 +0000 |
---|---|---|
committer | taw27 <taw27@bf6ca9ba-c028-0410-8290-897cf20841d1> | 2007-09-07 22:32:35 +0000 |
commit | ece9b68790622809b429bb54a124421e479c75c9 (patch) | |
tree | 64504947b5f81baa7866c37658b8a6fcf446916c /src | |
parent | cec71f84b2db92d6c5af4b2b28cbf294e2cc948f (diff) |
Add eFOM (needs debugging)
git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@127 bf6ca9ba-c028-0410-8290-897cf20841d1
Diffstat (limited to 'src')
-rw-r--r-- | src/ipr.c | 92 |
1 files changed, 83 insertions, 9 deletions
@@ -125,9 +125,7 @@ static void ipr_try_compact(Basis *basis) { } -/* Choose an initial set of three vectors for the expansion. - * This would probably be just as easy for the user to do... */ -static Basis *ipr_choose_initial_basis(ControlContext *ctx) { +static Basis *ipr_choose_random_basis(ControlContext *ctx) { Basis *basis; Reflection *reflection; @@ -135,7 +133,6 @@ static Basis *ipr_choose_initial_basis(ControlContext *ctx) { basis = malloc(sizeof(Basis)); - /* Determine a vaguely sensible starting basis */ basis->a.x = 0; basis->a.y = 0; basis->a.z = 0; basis->a.modulus = 1e30; basis->b.x = 0; basis->b.y = 0; basis->b.z = 0; basis->b.modulus = 1e30; basis->c.x = 0; basis->c.y = 0; basis->c.z = 0; basis->c.modulus = 1e30; @@ -191,7 +188,7 @@ static Basis *ipr_choose_initial_basis(ControlContext *ctx) { /* End of list? */ if ( !reflection ) { if ( (basis->a.modulus == 1e30) || (basis->b.modulus == 1e30) || (basis->c.modulus == 1e30) ) { - printf("IP: Failed to find appropriate starting basis!\n"); + //printf("IP: Failed to find appropriate starting basis!\n"); free(basis); return NULL; } @@ -201,17 +198,94 @@ static Basis *ipr_choose_initial_basis(ControlContext *ctx) { ipr_try_compact(basis); - basis->a.x = 1.842e9; basis->a.y = 0.0; basis->a.z = 0.0; - basis->b.x = 0.0; basis->b.y = 1.842e9; basis->b.z = 0.0; - basis->c.x = 0.0; basis->c.y = 0.0; basis->c.z = 1.842e9; +// basis->a.x = 1.842e9; basis->a.y = 0.0; basis->a.z = 0.0; +// basis->b.x = 0.0; basis->b.y = 1.842e9; basis->b.z = 0.0; +// basis->c.x = 0.0; basis->c.y = 0.0; basis->c.z = 1.842e9; + + return basis; + +} + +static double ipr_efom(ReflectionContext *rctx, Basis *basis) { + + int n_indexed, n_counted; + Reflection *cur; + + cur = rctx->reflections; + n_indexed = 0; + n_counted = 0; + while ( cur ) { + + if ( cur->type == REFLECTION_NORMAL ) { + + /* Can this basis approximately account for this reflection? */ + double n; + double p, q, r, u, v, w, d, e, f; + double h, k, l; + + /* Set up the coordinate transform from hkl to xyz */ + p = basis->a.x; q = basis->a.y; r = basis->a.z; + u = basis->b.x; v = basis->b.y; w = basis->b.z; + d = basis->c.x; e = basis->c.y; f = basis->c.z; + + /* Invert the matrix to get hkl from xyz */ + n = p*(v*f - w*e) - q*(u*f - w*d) + r*(u*e - d*v); + h = (p*(v*f-w*e)*cur->x + -q*(u*f-w*d)*cur->y + r*(u*e-v*d)*cur->z) / n; + k = (-u*(q*f-r*e)*cur->x + v*(p*f-r*d)*cur->y + -w*(p*e-q*d)*cur->z) / n; + l = (d*(q*w-r*v)*cur->x + -e*(p*w-r*u)*cur->y + f*(p*v-q*u)*cur->z) / n; + + /* Calculate the deviations in terms of |a|, |b| and |c| */ + h -= floor(h); k -= floor(k); l -= floor(l); + if ( h*h + k*k + l*l <= 0.1*0.1*0.1 ) n_indexed++; + + n_counted++; + + } + + cur = cur->next; + + } + + return (double)n_indexed / n_counted; +} + +static Basis *ipr_choose_initial_basis(ControlContext *ctx) { + + Basis *basis; + Basis *best_basis; + double best_efom; + int i; + + /* Track the best basis throughout the whole procedure */ + best_basis = NULL; best_efom = 0; + + printf("IP: Finding initial basis using eFOM...\n"); + for ( i=1; i<=1000; i++ ) { + + double efom; + + do { + basis = ipr_choose_random_basis(ctx); + } while ( !basis ); + efom = ipr_efom(ctx->reflectionctx, basis); + + if ( (efom > best_efom) || !best_basis ) { + if ( best_basis ) free(best_basis); + best_efom = efom; + best_basis = basis; + printf("IP: %3i: eFOM=%f\n", i, efom); + + } + + } + printf("IP: Initial basis:\n"); printf("IP: a = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->a.x, basis->a.y, basis->a.z, basis->a.modulus); printf("IP: b = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->b.x, basis->b.y, basis->b.z, basis->b.modulus); printf("IP: c = [%10.2e %10.2e %10.2e] mod=%10.2e\n", basis->c.x, basis->c.y, basis->c.z, basis->c.modulus); return basis; - } /* Generate list of reflections to analyse */ |