From b62612cc17eda59a3aa347d8a51e7ea235ad716b Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 19 Jun 2019 10:55:56 +0200 Subject: partialator: WIP on new minimiser --- src/partialator.c | 2 +- src/post-refinement.c | 610 +++++++++++++++----------------------------------- src/post-refinement.h | 2 +- 3 files changed, 183 insertions(+), 431 deletions(-) diff --git a/src/partialator.c b/src/partialator.c index f5d1a210..92fcfbe4 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -1500,7 +1500,7 @@ int main(int argc, char *argv[]) if ( !no_pr ) { refine_all(crystals, n_crystals, full, nthreads, pmodel, - 0, i+1, no_logs, sym, amb, scaleflags); + i+1, no_logs, sym, amb, scaleflags); } /* Create new reference if needed */ diff --git a/src/post-refinement.c b/src/post-refinement.c index 01ba8b14..b781d395 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -33,8 +33,6 @@ #include #include -#include -#include #include "image.h" #include "post-refinement.h" @@ -47,12 +45,27 @@ #include "scaling.h" #include "merge.h" +struct rf_alteration +{ + double rot_x; + double rot_y; + double delta_R; + double delta_wave; +}; + -struct prdata +struct rf_priv { - int refined; + const Crystal *cr; /**< Original crystal (before any refinement) */ + const RefList *full; + int serial; + int scaleflags; + + Crystal *cr_tgt; /**< Crystal to use for testing modifications */ + struct image image_tgt; /**< Image structure to go with cr_tgt */ }; + const char *str_prflag(enum prflag flag) { switch ( flag ) { @@ -84,16 +97,17 @@ const char *str_prflag(enum prflag flag) } -static void rotate_cell_xy(UnitCell *cell, double ang1, double ang2) +static void rotate_cell_xy(const UnitCell *source, UnitCell *tgt, + double ang1, double ang2) { double asx, asy, asz; double bsx, bsy, bsz; double csx, csy, csz; double xnew, ynew, znew; - cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); + cell_get_reciprocal(source, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); /* "a" around x */ xnew = asx; @@ -131,219 +145,55 @@ static void rotate_cell_xy(UnitCell *cell, double ang1, double ang2) znew = -csx*sin(ang2) + csz*cos(ang2); csx = xnew; csy = ynew; csz = znew; - cell_set_reciprocal(cell, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); -} - - -static const char *get_label(enum gparam p) -{ - switch ( p ) { - case GPARAM_ANG1 : return "First angle (radians)"; - case GPARAM_ANG2 : return "Second angle (radians)"; - case GPARAM_R : return "Profile radius (m^-1)"; - case GPARAM_WAVELENGTH : return "Wavelength (Angstroms)"; - default : return "unknown"; - } -} - - -/* We set all the step sizes to 1, then scale them. - * This way, the size of the simplex stays meaningful and we possibly also - * avoid some roundoff errors */ -static double get_scale(enum gparam p) -{ - switch ( p ) { - - case GPARAM_ANG1 : return deg2rad(0.05); - case GPARAM_ANG2 : return deg2rad(0.05); - case GPARAM_R : return 0.0005e9; - case GPARAM_WAVELENGTH : return 0.001e-10; - - default : return 0.0; - - } -} - - -struct rf_priv { - const Crystal *cr; - const RefList *full; - enum gparam rv[32]; - int verbose; - int serial; - gsl_vector *initial; - int scaleflags; - - /* For freeing later */ - gsl_vector *vals; - gsl_vector *step; - - /* So that it stays around until the end of minimisation */ - gsl_multimin_function f; -}; - - -static double get_actual_val(const gsl_vector *v, const gsl_vector *initial, - enum gparam *rv, int i) -{ - return gsl_vector_get(v, i) * get_scale(rv[i]) - + gsl_vector_get(initial, i); + cell_set_reciprocal(tgt, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); } -static void apply_parameters(const gsl_vector *v, const gsl_vector *initial, - enum gparam *rv, Crystal *cr) +static void apply_parameters(const Crystal *cr_orig, Crystal *cr_tgt, + struct rf_alteration alter) { - int i; - double ang1, ang2, R, lambda; + double R, lambda; - /* Default parameters if not used in refinement */ - ang1 = 0.0; - ang2 = 0.0; - R = crystal_get_profile_radius(cr); - lambda = crystal_get_image(cr)->lambda; + R = crystal_get_profile_radius(cr_orig); + lambda = crystal_get_image_const(cr_orig)->lambda; - for ( i=0; isize; i++ ) { - - double val; - - val = get_actual_val(v, initial, rv, i); - - switch ( rv[i] ) { - - case GPARAM_ANG1 : - ang1 = val; - break; - - case GPARAM_ANG2 : - ang2 = val; - break; - - case GPARAM_R : - R = val; - break; - - case GPARAM_WAVELENGTH : - lambda = val; - break; - - default : - ERROR("Don't understand parameter %i\n", rv[i]); - break; - - } - } - - rotate_cell_xy(crystal_get_cell(cr), ang1, ang2); - crystal_set_profile_radius(cr, R); - crystal_get_image(cr)->lambda = lambda; + rotate_cell_xy(crystal_get_cell_const(cr_orig), crystal_get_cell(cr_tgt), + alter.rot_x, alter.rot_y); + crystal_set_profile_radius(cr_tgt, R+alter.delta_R); + crystal_get_image(cr_tgt)->lambda = lambda+alter.delta_wave; } -static double residual_f(const gsl_vector *v, void *pp) +static double calc_residual(struct rf_priv *pv, struct rf_alteration alter, + int free) { - struct rf_priv *pv = pp; - RefList *list; - struct image im; - Crystal *cr; - double res; - int i; + apply_parameters(pv->cr, pv->cr_tgt, alter); - for ( i=0; isize; i++ ) { - if ( gsl_vector_get(v, i) > 100.0 ) return GSL_NAN; - } - - cr = crystal_copy(pv->cr); - im = *crystal_get_image(cr); - crystal_set_image(cr, &im); - crystal_set_cell(cr, cell_new_from_cell(crystal_get_cell(cr))); - - apply_parameters(v, pv->initial, pv->rv, cr); - - if ( fabs(crystal_get_profile_radius(cr)) > 5e9 ) { - cell_free(crystal_get_cell(cr)); - crystal_free(cr); - if ( pv->verbose ) STATUS("radius > 5e9\n"); - return GSL_NAN; + if ( fabs(crystal_get_profile_radius(pv->cr_tgt)) > 5e9 ) { + STATUS("radius > 5e9\n"); + return NAN; } /* Can happen with grid scans and certain --force-radius values */ - if ( fabs(crystal_get_profile_radius(cr)) < 0.0000001e9 ) { - cell_free(crystal_get_cell(cr)); - crystal_free(cr); - if ( pv->verbose ) STATUS("radius very small\n"); - return GSL_NAN; + if ( fabs(crystal_get_profile_radius(pv->cr_tgt)) < 0.0000001e9 ) { + STATUS("radius very small\n"); + return NAN; } - if ( im.lambda <= 0.0 ) { - cell_free(crystal_get_cell(cr)); - crystal_free(cr); - if ( pv->verbose ) STATUS("lambda < 0\n"); - return GSL_NAN; + if ( crystal_get_image(pv->cr_tgt)->lambda <= 0.0 ) { + STATUS("lambda < 0\n"); + return NAN; } - list = copy_reflist(crystal_get_reflections(cr)); - crystal_set_reflections(cr, list); - - update_predictions(cr); - calculate_partialities(cr, PMODEL_XSPHERE); + update_predictions(pv->cr_tgt); + calculate_partialities(pv->cr_tgt, PMODEL_XSPHERE); - if ( scale_one_crystal(cr, pv->full, pv->scaleflags) ) { - cell_free(crystal_get_cell(cr)); - reflist_free(crystal_get_reflections(cr)); - crystal_free(cr); - if ( pv->verbose ) STATUS("Bad scaling\n"); - return GSL_NAN; + if ( scale_one_crystal(pv->cr_tgt, pv->full, pv->scaleflags) ) { + //STATUS("Bad scaling\n"); + return NAN; } - res = residual(cr, pv->full, 0, NULL, NULL); - - cell_free(crystal_get_cell(cr)); - reflist_free(crystal_get_reflections(cr)); - crystal_free(cr); - - return res; -} - -static double get_initial_param(Crystal *cr, enum gparam p) -{ - switch ( p ) { - - case GPARAM_ANG1 : return 0.0; - case GPARAM_ANG2 : return 0.0; - case GPARAM_R : return crystal_get_profile_radius(cr); - case GPARAM_WAVELENGTH : return crystal_get_image(cr)->lambda; - - default: return 0.0; - - } -} - - -static int check_angle_shifts(gsl_vector *cur, const gsl_vector *initial, - enum gparam *rv, struct rf_priv *residual_f_priv) -{ - int i = 0; - double ang = 0.0; - - while ( rv[i] != GPARAM_EOL ) { - if ( (rv[i] == GPARAM_ANG1) || (rv[i] == GPARAM_ANG2) ) { - ang += fabs(get_actual_val(cur, initial, rv, i)); - } - rv++; - } - - if ( rad2deg(ang) > 5.0 ) { - ERROR("More than 5 degrees total rotation!\n"); - residual_f_priv->verbose = 1; - double res = residual_f(cur, residual_f_priv); - STATUS("residual after rotation = %e\n", res); - residual_f_priv->verbose = 2; - res = residual_f(initial, residual_f_priv); - STATUS("residual before rotation = %e\n", res); - return 1; - } - return 0; + return residual(pv->cr_tgt, pv->full, free, NULL, NULL); } @@ -587,90 +437,87 @@ void write_specgraph(Crystal *crystal, const RefList *full, } -static gsl_multimin_fminimizer *setup_minimiser(Crystal *cr, const RefList *full, - int verbose, int serial, - int scaleflags, - struct rf_priv *priv) +static void write_angle_grid(Crystal *cr, const RefList *full, + signed int cycle, int serial, int scaleflags) { - gsl_multimin_fminimizer *min; - int n_params = 0; - int i, r; - - /* The parameters to be refined */ - priv->rv[n_params++] = GPARAM_ANG1; - priv->rv[n_params++] = GPARAM_ANG2; - priv->rv[n_params++] = GPARAM_R; - priv->rv[n_params++] = GPARAM_WAVELENGTH; - priv->rv[n_params] = GPARAM_EOL; /* End of list */ - - priv->cr = cr; - priv->full = full; - priv->verbose = verbose; - priv->serial = serial; - priv->scaleflags = scaleflags; - - priv->f.f = residual_f; - priv->f.n = n_params; - priv->f.params = priv; - - priv->initial = gsl_vector_alloc(n_params); - priv->vals = gsl_vector_alloc(n_params); - priv->step = gsl_vector_alloc(n_params); - - for ( i=0; iinitial, i, get_initial_param(cr, priv->rv[i])); - gsl_vector_set(priv->vals, i, 0.0); - gsl_vector_set(priv->step, i, 1.0); - } + FILE *fh; + char fn[64]; + char ins[16]; + struct rf_priv priv; + RefList *list; + UnitCell *cell; - min = gsl_multimin_fminimizer_alloc(gsl_multimin_fminimizer_nmsimplex2, - n_params); - if ( min == NULL ) { - ERROR("Failed to allocate minimiser\n"); - gsl_vector_free(priv->vals); - gsl_vector_free(priv->step); - gsl_vector_free(priv->initial); - return NULL; + priv.cr = cr; + priv.full = full; + priv.serial = serial; + priv.scaleflags = scaleflags; + priv.cr_tgt = crystal_copy(cr); + priv.image_tgt = *crystal_get_image(cr); + crystal_set_image(priv.cr_tgt, &priv.image_tgt); + list = copy_reflist(crystal_get_reflections(cr)); + crystal_set_reflections(priv.cr_tgt, list); + cell = cell_new_from_cell(crystal_get_cell(cr)); + crystal_set_cell(priv.cr_tgt, cell); + + if ( cycle >= 0 ) { + snprintf(ins, 16, "%i", cycle); + } else { + ins[0] = 'F'; + ins[1] = '\0'; } - r = gsl_multimin_fminimizer_set(min, &priv->f, priv->vals, priv->step); - if ( r != 0 ) { - gsl_multimin_fminimizer_free(min); - gsl_vector_free(priv->vals); - gsl_vector_free(priv->step); - gsl_vector_free(priv->initial); - ERROR("Failed to set up minimiser: %s\n", gsl_strerror(r)); - return NULL; + snprintf(fn, 64, "pr-logs/grid-crystal%i-cycle%s-ang1-ang2.dat", + serial, ins); + fh = fopen(fn, "w"); + if ( fh != NULL ) { + double v1, v2; + fprintf(fh, "%e %e %e %s\n", -5.0e-3, 5.0e-3, 0.0, "ang1/rad"); + fprintf(fh, "%e %e %e %s\n", -5.0e-3, 5.0e-3, 0.0, "ang2/rad"); + for ( v2=-5.0e-3; v2<=5.0e-3; v2+=0.25e-3 ) { + int first=1; + for ( v1=-5.0e-3; v1<=5.0e-3; v1+=0.25e-3 ) { + double res; + struct rf_alteration alter; + alter.rot_x = v1; + alter.rot_y = v2; + alter.delta_R = 0.0; + alter.delta_wave = 0.0; + res = calc_residual(&priv, alter, 0); + if ( !first ) fprintf(fh, " "); + first = 0; + fprintf(fh, "%e", res); + } + fprintf(fh, "\n"); + } } + fclose(fh); - return min; + reflist_free(crystal_get_reflections(priv.cr_tgt)); + crystal_free(priv.cr_tgt); } -static void write_grid(Crystal *cr, const RefList *full, - signed int cycle, int serial, int scaleflags, - const enum gparam par1, const enum gparam par2, - const char *name) +static void write_radius_grid(Crystal *cr, const RefList *full, + signed int cycle, int serial, int scaleflags) { FILE *fh; char fn[64]; char ins[16]; - gsl_multimin_fminimizer *min; struct rf_priv priv; - int idx1, idx2; - int i; - - min = setup_minimiser(cr, full, 0, serial, scaleflags, &priv); - if ( min == NULL ) return; + RefList *list; + UnitCell *cell; - idx1 = 99; - idx2 = 99; - for ( i=0; i= 0 ) { snprintf(ins, 16, "%i", cycle); @@ -679,28 +526,23 @@ static void write_grid(Crystal *cr, const RefList *full, ins[1] = '\0'; } - snprintf(fn, 64, "pr-logs/grid-crystal%i-cycle%s-%s.dat", - serial, ins, name); + snprintf(fn, 64, "pr-logs/grid-crystal%i-cycle%s-R-wave.dat", + serial, ins); fh = fopen(fn, "w"); if ( fh != NULL ) { double v1, v2; - fprintf(fh, "%e %e %e %s\n", - -5.0*get_scale(par1)+get_initial_param(cr, par1), - 5.0*get_scale(par1)+get_initial_param(cr, par1), - get_initial_param(cr, par1), - get_label(par1)); - fprintf(fh, "%e %e %e %s\n", - -5.0*get_scale(par2)+get_initial_param(cr, par2), - 5.0*get_scale(par2)+get_initial_param(cr, par2), - get_initial_param(cr, par2), - get_label(par2)); - for ( v2=-5.0; v2<=5.0; v2+=0.25 ) { + fprintf(fh, "%e %e %e %s\n", -4e-13, 4e-13, 0.0, "wavelength change/m"); + fprintf(fh, "%e %e %e %s\n", -2e6, 2e6, 0.0, "radius change/m^-1"); + for ( v2=-2e6; v2<=2e6; v2+=40000 ) { int first=1; - for ( v1=-5.0; v1<=5.0; v1+=0.25 ) { + for ( v1=-4e-13; v1<=4e-13; v1+=8e-15 ) { double res; - gsl_vector_set(min->x, idx1, v1); - gsl_vector_set(min->x, idx2, v2); - res = residual_f(min->x, &priv); + struct rf_alteration alter; + alter.rot_x = 0.0; + alter.rot_y = 0.0; + alter.delta_R = v2; + alter.delta_wave = v1; + res = calc_residual(&priv, alter, 0); if ( !first ) fprintf(fh, " "); first = 0; fprintf(fh, "%e", res); @@ -710,36 +552,32 @@ static void write_grid(Crystal *cr, const RefList *full, } fclose(fh); - gsl_multimin_fminimizer_free(min); - gsl_vector_free(priv.initial); - gsl_vector_free(priv.vals); - gsl_vector_free(priv.step); + reflist_free(crystal_get_reflections(priv.cr_tgt)); + crystal_free(priv.cr_tgt); } void write_gridscan(Crystal *cr, const RefList *full, signed int cycle, int serial, int scaleflags) { - write_grid(cr, full, cycle, serial, scaleflags, - GPARAM_ANG1, GPARAM_ANG2, "ang1-ang2"); - write_grid(cr, full, cycle, serial, scaleflags, - GPARAM_ANG1, GPARAM_WAVELENGTH, "ang1-wave"); - write_grid(cr, full, cycle, serial, scaleflags, - GPARAM_R, GPARAM_WAVELENGTH, "R-wave"); + write_angle_grid(cr, full, cycle, serial, scaleflags); + write_radius_grid(cr, full, cycle, serial, scaleflags); } static void do_pr_refine(Crystal *cr, const RefList *full, - PartialityModel pmodel, int verbose, int serial, + PartialityModel pmodel, int serial, int cycle, int write_logs, SymOpList *sym, SymOpList *amb, int scaleflags) { - gsl_multimin_fminimizer *min; struct rf_priv priv; + struct rf_alteration alter; int n_iter = 0; - int status; - double residual_start, residual_free_start; + int status = 0; + double fom, freefom; + RefList *list; FILE *fh = NULL; + UnitCell *cell; try_reindex(cr, full, sym, amb, scaleflags); @@ -747,31 +585,26 @@ static void do_pr_refine(Crystal *cr, const RefList *full, ERROR("Bad scaling at start of refinement.\n"); return; } - residual_start = residual(cr, full, 0, NULL, NULL); - residual_free_start = residual(cr, full, 1, NULL, NULL); - if ( verbose ) { - STATUS("\nPR initial: dev = %10.5e, free dev = %10.5e\n", - residual_start, residual_free_start); - } + alter.rot_x = 0.0; + alter.rot_y = 0.0; + alter.delta_R = 0.0; + alter.delta_wave = 0.0; - min = setup_minimiser(cr, full, verbose, serial, scaleflags, &priv); - if ( min == NULL ) return; - - if ( verbose ) { - double res = residual_f(min->x, &priv); - double size = gsl_multimin_fminimizer_size(min); - STATUS("At start: %f %f %f %f ----> %f %f %e %f residual = %e size %f\n", - gsl_vector_get(min->x, 0), - gsl_vector_get(min->x, 1), - gsl_vector_get(min->x, 2), - gsl_vector_get(min->x, 3), - rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)), - rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)), - get_actual_val(min->x, priv.initial, priv.rv, 2), - get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10, - res, size); - } + priv.cr = cr; + priv.full = full; + priv.serial = serial; + priv.scaleflags = scaleflags; + priv.cr_tgt = crystal_copy(cr); + priv.image_tgt = *crystal_get_image(cr); + crystal_set_image(priv.cr_tgt, &priv.image_tgt); + list = copy_reflist(crystal_get_reflections(cr)); + crystal_set_reflections(priv.cr_tgt, list); + cell = cell_new_from_cell(crystal_get_cell(cr)); + crystal_set_cell(priv.cr_tgt, cell); + + fom = calc_residual(&priv, alter, 0); + freefom = calc_residual(&priv, alter, 1); if ( write_logs ) { @@ -780,116 +613,46 @@ static void do_pr_refine(Crystal *cr, const RefList *full, snprintf(fn, 63, "pr-logs/crystal%i-cycle%i.log", serial, cycle); fh = fopen(fn, "w"); if ( fh != NULL ) { - fprintf(fh, "iteration RtoReference CCtoReference nref " - "ang1 ang2 radius wavelength\n"); - double res = residual_f(min->x, &priv); - fprintf(fh, "%5i %10.8f %10.8f %5i %10.8f %10.8f %e %e\n", - n_iter, res, 0.0, 0, - rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)), - rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)), - get_actual_val(min->x, priv.initial, priv.rv, 2), - get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10); + fprintf(fh, "iter FoM FreeFoM ang1 " + "ang2 radius wavelength\n"); + fprintf(fh, "%5i %10.8f %10.8f %10.8f %10.8f %e %e\n", + n_iter, fom, freefom, + rad2deg(alter.rot_x), rad2deg(alter.rot_y), + crystal_get_profile_radius(cr)+alter.delta_R, + crystal_get_image(cr)->lambda+alter.delta_wave); } } do { - double res; - n_iter++; - status = gsl_multimin_fminimizer_iterate(min); - if ( status ) break; + /* FIXME: Actually minimise */ - res = residual_f(min->x, &priv); - if ( isnan(res) ) { - status = GSL_ENOPROG; - break; - } - - if ( verbose ) { - double res = residual_f(min->x, &priv); - double size = gsl_multimin_fminimizer_size(min); - STATUS("%f %f %f %f ----> %f %f %e %f residual = %e size %f\n", - gsl_vector_get(min->x, 0), - gsl_vector_get(min->x, 1), - gsl_vector_get(min->x, 2), - gsl_vector_get(min->x, 3), - rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)), - rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)), - get_actual_val(min->x, priv.initial, priv.rv, 2), - get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10, - res, size); - } + fom = calc_residual(&priv, alter, 0); + freefom = calc_residual(&priv, alter, 1); if ( fh != NULL ) { - fprintf(fh, "%5i %10.8f %10.8f %5i %10.8f %10.8f %e %e\n", - n_iter, res, 0.0, 0, - rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)), - rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)), - get_actual_val(min->x, priv.initial, priv.rv, 2), - get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10); + fprintf(fh, "%5i %10.8f %10.8f %10.8f %10.8f %e %e\n", + n_iter, fom, freefom, + rad2deg(alter.rot_x), rad2deg(alter.rot_y), + crystal_get_profile_radius(cr)+alter.delta_R, + crystal_get_image(cr)->lambda+alter.delta_wave); } - status = gsl_multimin_test_size(min->size, 0.005); - - } while ( status == GSL_CONTINUE && n_iter < 1000 ); + } while ( n_iter < 10 ); /* FIXME: Termination criteria */ - if ( verbose ) { - STATUS("Done with refinement after %i iter\n", n_iter); - STATUS("status = %i (%s)\n", status, gsl_strerror(status)); - } - - if ( status == GSL_SUCCESS ) { - - if ( check_angle_shifts(min->x, priv.initial, priv.rv, &priv) ) return; - - if ( verbose ) { - - double res = residual_f(min->x, &priv); - double size = gsl_multimin_fminimizer_size(min); - STATUS("At end: %f %f %f %f ----> %f %f %e %f residual = %e size %f\n", - gsl_vector_get(min->x, 0), - gsl_vector_get(min->x, 1), - gsl_vector_get(min->x, 2), - gsl_vector_get(min->x, 3), - rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)), - rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)), - get_actual_val(min->x, priv.initial, priv.rv, 2), - get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10, - res, size); - - } - - if ( fh != NULL ) { - double res = residual_f(min->x, &priv); - fprintf(fh, "%5i %10.8f %10.8f %5i %10.8f %10.8f %e %e\n", - n_iter, res, 0.0, 0, - rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)), - rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)), - get_actual_val(min->x, priv.initial, priv.rv, 2), - get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10); - } + if ( status == 0 ) { /* Apply the final shifts */ - apply_parameters(min->x, priv.initial, priv.rv, cr); + apply_parameters(cr, cr, alter); update_predictions(cr); calculate_partialities(cr, PMODEL_XSPHERE); scale_one_crystal(cr, full, scaleflags); - if ( verbose ) { - - STATUS("After applying final shifts:\n"); - STATUS("PR final: dev = %10.5e, free dev = %10.5e\n", - residual(cr, full, 0, NULL, NULL), - residual(cr, full, 1, NULL, NULL)); - STATUS("Final R = %e m^-1\n", crystal_get_profile_radius(cr)); - - } - } else { - ERROR("Bad refinement: crystal %i (%s) after %i iterations\n", - serial, gsl_strerror(status), n_iter); + ERROR("Bad refinement: crystal %i after %i iterations\n", + serial, n_iter); } if ( write_logs ) { @@ -899,17 +662,15 @@ static void do_pr_refine(Crystal *cr, const RefList *full, } if ( crystal_get_profile_radius(cr) > 5e9 ) { - ERROR("Very large radius: crystal %i\n", serial); + ERROR("WARNING: Very large radius: crystal %i\n", serial); } if ( fh != NULL ) { fclose(fh); } - gsl_multimin_fminimizer_free(min); - gsl_vector_free(priv.initial); - gsl_vector_free(priv.vals); - gsl_vector_free(priv.step); + reflist_free(crystal_get_reflections(priv.cr_tgt)); + crystal_free(priv.cr_tgt); } @@ -919,8 +680,6 @@ struct refine_args Crystal *crystal; PartialityModel pmodel; int serial; - struct prdata prdata; - int verbose; int cycle; int no_logs; SymOpList *sym; @@ -946,15 +705,10 @@ static void refine_image(void *task, int id) int write_logs = 0; write_logs = !pargs->no_logs && (pargs->serial % 20 == 0); - pargs->prdata.refined = 0; - do_pr_refine(cr, pargs->full, pargs->pmodel, pargs->verbose, + do_pr_refine(cr, pargs->full, pargs->pmodel, pargs->serial, pargs->cycle, write_logs, pargs->sym, pargs->amb, pargs->scaleflags); - - if ( crystal_get_user_flag(cr) == 0 ) { - pargs->prdata.refined = 1; - } } @@ -988,7 +742,7 @@ static void done_image(void *vqargs, void *task) void refine_all(Crystal **crystals, int n_crystals, RefList *full, int nthreads, PartialityModel pmodel, - int verbose, int cycle, int no_logs, + int cycle, int no_logs, SymOpList *sym, SymOpList *amb, int scaleflags) { struct refine_args task_defaults; @@ -997,8 +751,6 @@ void refine_all(Crystal **crystals, int n_crystals, task_defaults.full = full; task_defaults.crystal = NULL; task_defaults.pmodel = pmodel; - task_defaults.prdata.refined = 0; - task_defaults.verbose = verbose; task_defaults.cycle = cycle; task_defaults.no_logs = no_logs; task_defaults.sym = sym; diff --git a/src/post-refinement.h b/src/post-refinement.h index e2873b76..0c84f164 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -60,7 +60,7 @@ extern const char *str_prflag(enum prflag flag); extern void refine_all(Crystal **crystals, int n_crystals, RefList *full, int nthreads, PartialityModel pmodel, - int verbose, int cycle, int no_logs, + int cycle, int no_logs, SymOpList *sym, SymOpList *amb, int scaleflags); extern void write_gridscan(Crystal *cr, const RefList *full, -- cgit v1.2.3