diff options
Diffstat (limited to 'src')
-rw-r--r-- | src/refine.c | 76 | ||||
-rw-r--r-- | src/refine.h | 2 |
2 files changed, 54 insertions, 24 deletions
diff --git a/src/refine.c b/src/refine.c index 4b6793e..4d2a12e 100644 --- a/src/refine.c +++ b/src/refine.c @@ -25,6 +25,7 @@ #include "image.h" #include "reproject.h" #include "mapping.h" +#include "refine.h" /* A simplex is an array of ten of these */ typedef struct { @@ -40,7 +41,26 @@ typedef struct { void refine_do_sequence(ControlContext *ctx) { - + double omega_offs; + + for ( omega_offs=-2.0; omega_offs<=2.0; omega_offs+=0.01 ) { + + double fit; + int i; + + for ( i=0; i<ctx->images->n_images; i++ ) { + ctx->images->images[i].omega += omega_offs; + } + + fit = refine_do_cell(ctx); + printf("RF: omega_offs=%f, fit=%f nm^-1\n", omega_offs, fit/1e9); + + for ( i=0; i<ctx->images->n_images; i++ ) { + ctx->images->images[i].omega -= omega_offs; + } + + } + } static double refine_mean_dev(Deviation *d, int nf, SimplexVertex *s, int i) { @@ -67,6 +87,7 @@ static double refine_mean_dev(Deviation *d, int nf, SimplexVertex *s, int i) { } +#if 0 static void refine_display_simplex(SimplexVertex *s) { int i; @@ -79,6 +100,7 @@ static void refine_display_simplex(SimplexVertex *s) { } } +#endif /* Expand the simplex across from vertex v_worst by factor 'fac'. * fac = -1 is a reflection @@ -103,14 +125,13 @@ static void refine_simplex_transform(SimplexVertex *s, int v_worst, double fac) } -static void refine_iteration(SimplexVertex *s, Deviation *d, int nf) { +static double refine_iteration(SimplexVertex *s, Deviation *d, int nf) { int v_worst, v_best, v_second_worst, i; double fom_worst, fom_new, fom_best, fom_second_worst; - /* Start an iteration */ - refine_display_simplex(s); - +// printf("---\n"); + /* Find the least favourable vertex of the simplex */ v_worst = 0; fom_worst = 0.0; @@ -124,43 +145,40 @@ static void refine_iteration(SimplexVertex *s, Deviation *d, int nf) { fom = refine_mean_dev(d, nf, s, i); - printf("Mean deviation at simplex vertex %i = %f nm^-1\n", i, fom/1e9); + //printf("Mean deviation at simplex vertex %i = %f nm^-1\n", i, fom/1e9); if ( fom > fom_worst ) { v_second_worst = v_worst; fom_second_worst = fom_worst; fom_worst = fom; v_worst = i; } - if ( fom > fom_best ) { + if ( fom < fom_best ) { fom_best = fom; v_best = i; } } - printf("The worst vertex is number %i\n", v_worst); - printf("The second worst vertex is number %i\n", v_second_worst); - printf("The best vertex is number %i\n", v_best); /* Reflect this vertex across the opposite face */ refine_simplex_transform(s, v_worst, -1.0); - + /* Is the worst vertex any better? */ fom_new = refine_mean_dev(d, nf, s, v_worst); - printf("New mean deviation for the worst vertex after reflection is %f nm^-1\n", fom_new/1e9); +// printf("New mean deviation for the worst vertex after reflection is %f nm^-1\n", fom_new/1e9); if ( fom_new > fom_worst ) { double fom_new_new; /* It's worse than before. Contract in 1D and see if that helps. */ - printf("Worse. Trying a 1D contraction...\n"); +// printf("Worse. Trying a 1D contraction...\n"); refine_simplex_transform(s, v_worst, 0.5); fom_new_new = refine_mean_dev(d, nf, s, v_worst); - printf("Mean deviation after 1D contraction is %f nm^-1\n", fom_new_new/1e9); +// printf("Mean deviation after 1D contraction is %f nm^-1\n", fom_new_new/1e9); if ( fom_new_new > fom_second_worst ) { int i; - printf("Not as good as the second worst vertex: contracting around the best vertex (%i)\n", v_best); +// printf("Not as good as the second worst vertex: contracting around the best vertex (%i)\n", v_best); for ( i=0; i<10; i++ ) { if ( i != v_best ) refine_simplex_transform(s, i, 0.5); } @@ -173,38 +191,41 @@ static void refine_iteration(SimplexVertex *s, Deviation *d, int nf) { double fom_new_new; SimplexVertex save; - printf("This is better. Trying to expand...\n"); +// printf("This is better. Trying to expand...\n"); save = s[v_worst]; refine_simplex_transform(s, v_worst, 2); /* Better? */ fom_new_new = refine_mean_dev(d, nf, s, v_worst); - printf("Mean deviation after expansion is %f nm^-1\n", fom_new_new/1e9); +// printf("Mean deviation after expansion is %f nm^-1\n", fom_new_new/1e9); if ( fom_new_new > fom_new ) { /* "Got too confident" */ s[v_worst] = save; - printf("Got too confident - reverting\n"); +// printf("Got too confident - reverting\n"); } /* else yay. */ } + + return fom_worst - fom_best; } -void refine_do_cell(ControlContext *ctx) { +double refine_do_cell(ControlContext *ctx) { SimplexVertex s[10]; Deviation *d; double delta; int i, nf, f, it, maxiter; + const double tol = 0.001e9; /* Stopping condition */ if ( !ctx->cell_lattice ) { displaywindow_error("No reciprocal unit cell has been found.", ctx->dw); - return; + return -1; } if ( ctx->images->n_images == 0 ) { displaywindow_error("There are no images to refine against.", ctx->dw); - return; + return -1; } /* Determine the size of the 'deviation table' */ @@ -291,13 +312,22 @@ void refine_do_cell(ControlContext *ctx) { /* Iterate */ maxiter = 500; for ( it=0; it<maxiter; it++ ) { - printf("Simplex method iteration %i\n", it); - refine_iteration(s, d, nf); + + double conv; + +// printf("Simplex method iteration %i\n", it); + conv = refine_iteration(s, d, nf); + if ( conv < tol ) { + printf("Converged after %i iterations (%f)\n", it, conv); + break; + } } ctx->images->images[ctx->dw->cur_image].rflist = NULL; reproject_lattice_changed(ctx); displaywindow_update(ctx->dw); + + return refine_mean_dev(d, nf, s, 0); } diff --git a/src/refine.h b/src/refine.h index d3ca6d9..38dee22 100644 --- a/src/refine.h +++ b/src/refine.h @@ -19,7 +19,7 @@ #include "control.h" extern void refine_do_sequence(ControlContext *ctx); -extern void refine_do_cell(ControlContext *ctx); +extern double refine_do_cell(ControlContext *ctx); #endif /* REFINE_H */ |