From a0aae9ef83b3c8c2b223f933fb11e9774dd440fc Mon Sep 17 00:00:00 2001 From: taw27 Date: Fri, 25 Apr 2008 16:06:33 +0000 Subject: Tidy up refinement unit test git-svn-id: svn://cook.msm.cam.ac.uk:745/diff-tomo/dtr@276 bf6ca9ba-c028-0410-8290-897cf20841d1 --- src/refine.c | 65 ++++++++++++++++++++++++++++++++++++++++-------------------- 1 file changed, 43 insertions(+), 22 deletions(-) (limited to 'src/refine.c') diff --git a/src/refine.c b/src/refine.c index e92939d..d0fc8cc 100644 --- a/src/refine.c +++ b/src/refine.c @@ -29,8 +29,15 @@ #include "gtk-valuegraph.h" #include "utils.h" +/* Divide numbers by this for display */ #define DISPFACTOR 1.0e9 +/* Number of parameters */ +#define NUM_PARAMS 9 + +/* Refine debug */ +#define REFINE_DEBUG 0 + /* A simplex is an array of ten of these */ typedef struct { double dax; double dbx; double dcx; @@ -157,7 +164,7 @@ static void refine_simplex_transform(SimplexVertex *s, int v_worst, double fac) centre.dcx = 0.0; centre.dcy = 0.0; centre.dcz = 0.0; centre.does_nothing = 0; nv = 0; - for ( i=0; i<5; i++ ) { + for ( i=0; i<=NUM_PARAMS; i++ ) { if ( i != v_worst ) { 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; @@ -186,8 +193,10 @@ 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]); + if ( REFINE_DEBUG ) { + printf(" After transformation: "); + refine_display_simplex(s[v_worst]); + } } @@ -206,7 +215,7 @@ static void refine_simplex_contract(SimplexVertex *s, int v, int v_best) { } -static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug) { +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; @@ -218,14 +227,14 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug fom_best = 100e9; v_second_worst = 0; fom_second_worst = 0.0; - if ( debug ) printf("Vertex FoM/nm^-1 DoesNothing\n"); - for ( i=0; i<5; i++ ) { + if ( REFINE_DEBUG ) printf("Vertex FoM/nm^-1 DoesNothing\n"); + for ( i=0; i<=NUM_PARAMS; i++ ) { double fom; fom = refine_mean_dev(d, nf, s, i); - if ( debug ) printf("%6i %8f %s\n", i, fom/DISPFACTOR, s[i].does_nothing?"*":" "); + if ( REFINE_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; @@ -238,29 +247,29 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug } } - if ( debug ) printf("The worst vertex is number %i\n", v_worst); + if ( REFINE_DEBUG ) printf("The worst vertex is number %i\n", v_worst); /* 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); - if ( debug ) printf("New mean deviation for the worst vertex after reflection is %f nm^-1\n", fom_new/DISPFACTOR); + if ( REFINE_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; /* It's worse than before. Contract in 1D and see if that helps. */ - if ( debug ) printf("Worse. Trying a 1D contraction...\n"); + if ( REFINE_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/DISPFACTOR); + if ( REFINE_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; - if ( debug ) printf("Not as good as the second worst vertex: contracting around the best vertex (%i)\n", v_best); - for ( i=0; i<5; i++ ) { + if ( REFINE_DEBUG ) printf("Not as good as the second worst vertex: contracting around the best vertex (%i)\n", v_best); + for ( i=0; i<=NUM_PARAMS; i++ ) { if ( i != v_best ) refine_simplex_contract(s, i, v_best); } @@ -272,19 +281,19 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug double fom_new_new; SimplexVertex save; - if ( debug ) printf("This is better. Trying to expand...\n"); + if ( REFINE_DEBUG ) printf("This is better. Trying to expand...\n"); save = s[v_worst]; 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/DISPFACTOR); + if ( REFINE_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; - if ( debug ) printf("Got too confident - reverting\n"); + if ( REFINE_DEBUG ) printf("Got too confident - reverting\n"); } else { - if ( debug ) printf("Better still. Great.\n"); + if ( REFINE_DEBUG ) printf("Better still. Great.\n"); } } else { @@ -294,6 +303,19 @@ static double refine_iteration(SimplexVertex *s, Deviation *d, int nf, int debug } + /* Check convergence and return */ + fom_worst = 0; + fom_best = 100e9; + for ( i=0; i<=NUM_PARAMS; i++ ) { + double fom; + fom = refine_mean_dev(d, nf, s, i); + if ( (s[i].does_nothing == 0) && (fom > fom_worst) ) { + fom_worst = fom; + } + if ( fom < fom_best ) { + fom_best = fom; + } + } return fom_worst - fom_best; } @@ -306,7 +328,6 @@ double refine_do_cell(ControlContext *ctx) { int i, nf, f, it, maxiter; const double tol = 0.0001e9; /* Stopping condition */ //const double tol = 0.001; /* For testing */ - int debug = 1; if ( !ctx->cell_lattice ) { displaywindow_error("No reciprocal unit cell has been found.", ctx->dw); @@ -332,7 +353,7 @@ double refine_do_cell(ControlContext *ctx) { } } - if ( debug ) printf("RF: There are %i partnered features in total\n", nf); + if ( REFINE_DEBUG ) printf("RF: There are %i partnered features in total\n", nf); /* Initialise the 'deviation table' */ d = malloc(nf*sizeof(Deviation)); @@ -410,10 +431,10 @@ double refine_do_cell(ControlContext *ctx) { // refine_display_simplex(s[i]); //} - if ( debug ) printf("------------------- Simplex method iteration %i -------------------\n", it); - conv = refine_iteration(s, d, nf, debug); + if ( REFINE_DEBUG ) printf("------------------- Simplex method iteration %i -------------------\n", it); + conv = refine_iteration(s, d, nf); if ( conv < tol ) { - if ( debug ) printf("RF: Converged after %i iterations (%f nm^-1)\n", it, conv/DISPFACTOR); + if ( REFINE_DEBUG ) printf("RF: Converged after %i iterations (%f nm^-1)\n", it, conv/DISPFACTOR); break; } -- cgit v1.2.3