From 5b828c42c58e48ac98d2c606df84a18cfd2652ae Mon Sep 17 00:00:00 2001 From: taw27 Date: Fri, 25 Apr 2008 11:57:43 +0000 Subject: Test minimisation works (yay!) git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@274 bf6ca9ba-c028-0410-8290-897cf20841d1 --- src/refine.c | 58 +++++++++++++++++++++++++++++++++------------------------- 1 file changed, 33 insertions(+), 25 deletions(-) diff --git a/src/refine.c b/src/refine.c index 45402d1..fc983ab 100644 --- a/src/refine.c +++ b/src/refine.c @@ -29,6 +29,8 @@ #include "gtk-valuegraph.h" #include "utils.h" +#define DISPFACTOR 1.0 + /* A simplex is an array of ten of these */ typedef struct { double dax; double dbx; double dcx; @@ -71,7 +73,7 @@ void refine_do_sequence(ControlContext *ctx) { reproject_lattice_changed(ctx); fit = refine_do_cell(ctx); - printf("RF: omega_offs=%f deg, fit=%f nm^-1\n", rad2deg(omega_offs), fit/1e9); + printf("RF: omega_offs=%f deg, fit=%f nm^-1\n", rad2deg(omega_offs), fit/DISPFACTOR); fit_vals[idx++] = fit; if ( fit < fit_best ) { fit_best = fit; @@ -94,7 +96,7 @@ void refine_do_sequence(ControlContext *ctx) { gtk_widget_show_all(window_fit); /* Perform final refinement */ - printf("Best omega offset = %f deg (%f nm^-1)\n", rad2deg(omega_offs_best), fit_best/1e9); + printf("Best omega offset = %f deg (%f nm^-1)\n", rad2deg(omega_offs_best), fit_best/DISPFACTOR); for ( j=0; jimages->n_images; j++ ) { ctx->images->images[j].omega += omega_offs_best; } @@ -124,18 +126,18 @@ static double refine_mean_dev(Deviation *d, int nf, SimplexVertex *s, int i) { } -// return fabs(s[i].dax-10.0) + fabs(s[i].day-20.0) + fabs(s[i].dbx-30.0) + fabs(s[i].dby-40.0); + return fabs(s[i].dax-10.0) + fabs(s[i].day-20.0) + fabs(s[i].dbx-30.0) + fabs(s[i].dby-40.0); - return fom/nf; +// return fom/nf; } static void refine_display_simplex(SimplexVertex sv) { - printf("%8f %8f %8f %8f %8f %8f %8f %8f %8f\n", - sv.dax/1e9, sv.day/1e9, sv.daz/1e9, - sv.dbx/1e9, sv.dby/1e9, sv.dbz/1e9, - sv.dcx/1e9, sv.dcy/1e9, sv.dcz/1e9); + printf("%9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f %9.6f\n", + sv.dax/DISPFACTOR, sv.day/DISPFACTOR, sv.daz/DISPFACTOR, + sv.dbx/DISPFACTOR, sv.dby/DISPFACTOR, sv.dbz/DISPFACTOR, + sv.dcx/DISPFACTOR, sv.dcy/DISPFACTOR, sv.dcz/DISPFACTOR); } @@ -153,12 +155,13 @@ static void refine_simplex_transform(SimplexVertex *s, int v_worst, double fac) centre.dax = 0.0; centre.day = 0.0; centre.daz = 0.0; centre.dbx = 0.0; centre.dby = 0.0; centre.dbz = 0.0; centre.dcx = 0.0; centre.dcy = 0.0; centre.dcz = 0.0; + centre.does_nothing = 0; nv = 0; for ( i=0; i<5; i++ ) { if ( i != v_worst ) { - centre.dax += s[i].dax; centre.dax += s[i].day; centre.dax += s[i].daz; - centre.dax += s[i].dbx; centre.dax += s[i].dby; centre.dax += s[i].dbz; - centre.dax += s[i].dcx; centre.dax += s[i].dcy; centre.dax += s[i].dcz; + centre.dax += s[i].dax; centre.day += s[i].day; centre.daz += s[i].daz; + centre.dbx += s[i].dbx; centre.dby += s[i].dby; centre.dbz += s[i].dbz; + centre.dcx += s[i].dcx; centre.dcy += s[i].dcy; centre.dcz += s[i].dcz; nv++; } } @@ -166,11 +169,11 @@ static void refine_simplex_transform(SimplexVertex *s, int v_worst, double fac) centre.dbx /= nv; centre.dby /= nv; centre.dbz /= nv; centre.dcx /= nv; centre.dcy /= nv; centre.dcz /= nv; -// printf(" Midpoint: "); -// refine_display_simplex(centre); + printf("Before transformation: "); + refine_display_simplex(s[v_worst]); -// printf("Before transformation: "); -// refine_display_simplex(s[v_worst]); + printf(" Midpoint: "); + refine_display_simplex(centre); /* Do the transformation */ s[v_worst].dax = centre.dax + fac * (s[v_worst].dax - centre.dax); @@ -183,8 +186,8 @@ static void refine_simplex_transform(SimplexVertex *s, int v_worst, double fac) s[v_worst].dcy = centre.dcy + fac * (s[v_worst].dcy - centre.dcy); s[v_worst].dcz = centre.dcz + fac * (s[v_worst].dcz - centre.dcz); -// printf(" After transformation: "); -// refine_display_simplex(s[v_worst]); + printf(" After transformation: "); + refine_display_simplex(s[v_worst]); } @@ -222,7 +225,7 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug fom = refine_mean_dev(d, nf, s, i); - if ( debug ) printf("%6i %8f %s\n", i, fom/1e9, s[i].does_nothing?"*":" "); + if ( debug ) printf("%6i %8f %s\n", i, fom/DISPFACTOR, s[i].does_nothing?"*":" "); if ( (s[i].does_nothing == 0) && (fom > fom_worst) ) { v_second_worst = v_worst; fom_second_worst = fom_worst; @@ -242,7 +245,7 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug /* Is the worst vertex any better? */ fom_new = refine_mean_dev(d, nf, s, v_worst); - if ( debug ) printf("New mean deviation for the worst vertex after reflection is %f nm^-1\n", fom_new/1e9); + if ( debug ) printf("New mean deviation for the worst vertex after reflection is %f nm^-1\n", fom_new/DISPFACTOR); if ( fom_new > fom_worst ) { double fom_new_new; @@ -251,7 +254,7 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug if ( debug ) printf("Worse. Trying a 1D contraction...\n"); refine_simplex_transform(s, v_worst, -0.5); /* Minus puts it back on the original side of the 'good' face */ fom_new_new = refine_mean_dev(d, nf, s, v_worst); - if ( debug ) printf("Mean deviation after 1D contraction is %f nm^-1\n", fom_new_new/1e9); + if ( debug ) printf("Mean deviation after 1D contraction is %f nm^-1\n", fom_new_new/DISPFACTOR); if ( fom_new_new > fom_second_worst ) { int i; @@ -275,7 +278,7 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug refine_simplex_transform(s, v_worst, 2.0); /* +ve means stay on this side of the 'good' face */ /* Better? */ fom_new_new = refine_mean_dev(d, nf, s, v_worst); - if ( debug ) printf("Mean deviation after expansion is %f nm^-1\n", fom_new_new/1e9); + if ( debug ) printf("Mean deviation after expansion is %f nm^-1\n", fom_new_new/DISPFACTOR); if ( fom_new_new > fom_new ) { /* "Got too confident" */ s[v_worst] = save; @@ -301,7 +304,8 @@ double refine_do_cell(ControlContext *ctx) { Deviation *d; double delta; int i, nf, f, it, maxiter; - const double tol = 0.001e9; /* Stopping condition */ + //const double tol = 0.001e9; /* Stopping condition */ + const double tol = 0.001; /* For testing */ int debug = 1; if ( !ctx->cell_lattice ) { @@ -367,7 +371,7 @@ double refine_do_cell(ControlContext *ctx) { mapping_rotate(dx, dy, 0.0, &dlx, &dly, &dlz, image->omega, image->tilt); rf->partner->x = old_x; rf->partner->y = old_y; - //printf("=> %+10.5f %+10.5f %+10.5f nm^-1\n", dlx/1e9, dly/1e9, dlz/1e9); + //printf("=> %+10.5f %+10.5f %+10.5f nm^-1\n", dlx/DISPFACTOR, dly/DISPFACTOR, dlz/DISPFACTOR); d[f].dx = dlx; d[f].dy = dly; @@ -381,7 +385,7 @@ double refine_do_cell(ControlContext *ctx) { assert( f == nf ); /* Initialise the simplex */ - delta = 10.0;//0.01e9; + delta = 5.0;//0.01e9; s[0].dax = 0.0; s[0].dbx = 0.0; s[0].dcx = 0.0; 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; @@ -402,10 +406,14 @@ double refine_do_cell(ControlContext *ctx) { double conv; + for ( i=0; i<5; i++ ) { + refine_display_simplex(s[i]); + } + if ( debug ) printf("------------------- Simplex method iteration %i -------------------\n", it); conv = refine_iteration(s, d, nf, debug); if ( conv < tol ) { - if ( debug ) printf("RF: Converged after %i iterations (%f nm^-1)\n", it, conv/1e9); + if ( debug ) printf("RF: Converged after %i iterations (%f nm^-1)\n", it, conv/DISPFACTOR); break; } -- cgit v1.2.3