diff options
-rw-r--r-- | src/refine.c | 150 |
1 files changed, 120 insertions, 30 deletions
diff --git a/src/refine.c b/src/refine.c index 7871dc1..6eff91d 100644 --- a/src/refine.c +++ b/src/refine.c @@ -63,11 +63,115 @@ static double refine_mean_dev(Deviation *d, int nf, SimplexVertex *s, int i) { } - return fom; + return fom/nf; } -static void refine_simplex_reflect(SimplexVertex *s, int v_worst) { +static void refine_display_simplex(SimplexVertex *s) { + + int i; + + for ( i=0; i<10; i++ ) { + printf("Vertex %i: %8f %8f %8f %8f %8f %8f %8f %8f %8f\n", + i, s[i].dax/1e9, s[i].day/1e9, s[i].daz/1e9, + s[i].dbx/1e9, s[i].dby/1e9, s[i].dbz/1e9, + s[i].dcx/1e9, s[i].dcy/1e9, s[i].dcz/1e9); + } + +} + +/* Expand the simplex across from vertex v_worst by factor 'fac'. + * fac = -1 is a reflection + * fac = +n is a 1d expansion + */ +static void refine_simplex_transform(SimplexVertex *s, int v_worst, double fac) { + + int v_face; /* A simplex vertex which is on the face being reflected across */ + + v_face = 0; + if ( v_worst == 0 ) v_face = 1; + + s[v_worst].dax = s[v_face].dax - fac * (s[v_worst].dax - s[v_face].dax); + s[v_worst].day = s[v_face].day - fac * (s[v_worst].day - s[v_face].dax); + s[v_worst].daz = s[v_face].daz - fac * (s[v_worst].daz - s[v_face].dax); + s[v_worst].dbx = s[v_face].dbx - fac * (s[v_worst].dbx - s[v_face].dbx); + s[v_worst].dby = s[v_face].dby - fac * (s[v_worst].dby - s[v_face].dbx); + s[v_worst].dbz = s[v_face].dbz - fac * (s[v_worst].dbz - s[v_face].dbx); + s[v_worst].dcx = s[v_face].dcx - fac * (s[v_worst].dcx - s[v_face].dcx); + s[v_worst].dcy = s[v_face].dcy - fac * (s[v_worst].dcy - s[v_face].dcx); + s[v_worst].dcz = s[v_face].dcz - fac * (s[v_worst].dcz - s[v_face].dcx); + +} + +static void refine_iteration(SimplexVertex *s, Deviation *d, int nf) { + + int v_worst, v_best, i; + double fom_worst, fom_new, fom_best; + + /* Start an iteration */ + refine_display_simplex(s); + + /* Find the least favourable vertex of the simplex */ + v_worst = 0; + fom_worst = 0.0; + v_best = 0; + fom_best = 100e9; + for ( i=0; i<10; i++ ) { + + double fom; + + fom = refine_mean_dev(d, nf, s, i); + + printf("Mean deviation at simplex vertex %i = %f nm^-1\n", i, fom/1e9); + if ( fom > fom_worst ) { + fom_worst = fom; + v_worst = i; + } + if ( fom > fom_best ) { + fom_best = fom; + v_best = i; + } + + } + printf("The worst vertex is number %i\n", v_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); + + /* 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); + if ( fom_new > fom_worst ) { + + /* It's worse than before. Contract around the best vertex */ + int i; + printf("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); + } + + } else { + + /* It's better. Try to expand in this direction */ + double fom_new_new; + SimplexVertex save; + + 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); + if ( fom_new_new > fom_new ) { + /* "Got too confident" */ + s[v_worst] = save; + printf("Got too confident - reverting\n"); + } /* else yay. */ + + } + } void refine_do_cell(ControlContext *ctx) { @@ -75,8 +179,7 @@ void refine_do_cell(ControlContext *ctx) { SimplexVertex s[10]; Deviation *d; double delta; - int i, nf, f, v_worst; - double fom_worst; + int i, nf, f, it, maxiter; if ( !ctx->cell_lattice ) { displaywindow_error("No reciprocal unit cell has been found.", ctx->dw); @@ -88,7 +191,7 @@ void refine_do_cell(ControlContext *ctx) { return; } - /* Create the table of indicies and deviations */ + /* Determine the size of the 'deviation table' */ nf = 0; for ( i=0; i<ctx->images->n_images; i++ ) { int j; @@ -104,6 +207,7 @@ void refine_do_cell(ControlContext *ctx) { } printf("RF: There are %i partnered features in total\n", nf); + /* Initialise the 'deviation table' */ d = malloc(nf*sizeof(Deviation)); f = 0; for ( i=0; i<ctx->images->n_images; i++ ) { @@ -159,35 +263,21 @@ void refine_do_cell(ControlContext *ctx) { s[0].day = 0.0; s[0].dby = 0.0; s[0].dcy = 0.0; s[0].daz = 0.0; s[0].dbz = 0.0; s[0].dcz = 0.0; memcpy(&s[1], &s[0], sizeof(SimplexVertex)); s[1].dax = delta; - memcpy(&s[2], &s[0], sizeof(SimplexVertex)); s[2].dbx = delta; - memcpy(&s[3], &s[0], sizeof(SimplexVertex)); s[3].dcx = delta; - memcpy(&s[4], &s[0], sizeof(SimplexVertex)); s[4].day = delta; + memcpy(&s[2], &s[0], sizeof(SimplexVertex)); s[2].day = delta; + memcpy(&s[3], &s[0], sizeof(SimplexVertex)); s[3].daz = delta; + memcpy(&s[4], &s[0], sizeof(SimplexVertex)); s[4].dbx = delta; memcpy(&s[5], &s[0], sizeof(SimplexVertex)); s[5].dby = delta; - memcpy(&s[6], &s[0], sizeof(SimplexVertex)); s[6].dcy = delta; - memcpy(&s[7], &s[0], sizeof(SimplexVertex)); s[7].daz = delta; - memcpy(&s[8], &s[0], sizeof(SimplexVertex)); s[8].dbz = delta; + memcpy(&s[6], &s[0], sizeof(SimplexVertex)); s[6].dbz = delta; + memcpy(&s[7], &s[0], sizeof(SimplexVertex)); s[7].dcx = delta; + memcpy(&s[8], &s[0], sizeof(SimplexVertex)); s[8].dcy = delta; memcpy(&s[9], &s[0], sizeof(SimplexVertex)); s[9].dcz = delta; - /* Find the least favourable vertex of the simplex */ - v_worst = 0; - fom_worst = 0.0; - for ( i=0; i<10; i++ ) { - - double fom; - - fom = refine_mean_dev(d, nf, s, i); - - printf("Mean deviation at simplex vertex %i = %f nm^-1\n", i, (fom/nf)/1e9); - if ( fom > fom_worst ) { - fom_worst = fom; - v_worst = i; - } - + /* Iterate */ + maxiter = 2; + for ( it=0; it<maxiter; it++ ) { + printf("Simplex method iteration %i\n", it); + refine_iteration(s, d, nf); } - printf("The worst vertex is number %i\n", v_worst); - - /* Reflect this vertex across the opposite face */ - refine_simplex_reflect(s, v_worst); ctx->images->images[ctx->dw->cur_image].rflist = NULL; reproject_lattice_changed(ctx); |