aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/ambigator.c39
-rw-r--r--src/geoptimiser.c384
-rw-r--r--src/im-sandbox.c14
-rw-r--r--src/indexamajig.c35
-rw-r--r--src/merge.c128
-rw-r--r--src/merge.h9
-rw-r--r--src/partial_sim.c5
-rw-r--r--src/partialator.c460
-rw-r--r--src/post-refinement.c1146
-rw-r--r--src/post-refinement.h20
-rw-r--r--src/process_hkl.c8
-rw-r--r--src/process_image.c2
-rw-r--r--src/rejection.c2
-rw-r--r--src/scaling.c461
-rw-r--r--src/scaling.h16
15 files changed, 1769 insertions, 960 deletions
diff --git a/src/ambigator.c b/src/ambigator.c
index ba7629ab..bbf1d06e 100644
--- a/src/ambigator.c
+++ b/src/ambigator.c
@@ -674,7 +674,8 @@ static void reindex_reflections(FILE *fh, FILE *ofh, int assignment,
/* This is nasty, but means the output includes absolutely everything in the
* input, even stuff ignored by read_chunk() */
static void write_reindexed_stream(const char *infile, const char *outfile,
- int *assignments, SymOpList *amb)
+ int *assignments, SymOpList *amb,
+ int argc, char *argv[])
{
FILE *fh;
FILE *ofh;
@@ -683,6 +684,7 @@ static void write_reindexed_stream(const char *infile, const char *outfile,
int have_as = 0;
int have_bs = 0;
int have_cs = 0;
+ int done = 0;
fh = fopen(infile, "r");
if ( fh == NULL ) {
@@ -696,6 +698,38 @@ static void write_reindexed_stream(const char *infile, const char *outfile,
return;
}
+ /* Copy the header */
+ do {
+
+ char line[1024];
+ char *rval;
+
+ rval = fgets(line, 1023, fh);
+ if ( rval == NULL ) {
+ ERROR("Failed to read stream audit info.\n");
+ return;
+ }
+
+ if ( strncmp(line, "-----", 5) == 0 ) {
+
+ done = 1;
+
+ /* Add our own header */
+ fprintf(ofh, "Re-indexed by ambigator "CRYSTFEL_VERSIONSTRING"\n");
+ if ( argc > 0 ) {
+ for ( i=0; i<argc; i++ ) {
+ if ( i > 0 ) fprintf(ofh, " ");
+ fprintf(ofh, "%s", argv[i]);
+ }
+ fprintf(ofh, "\n");
+ }
+
+ }
+
+ fputs(line, ofh);
+
+ } while ( !done );
+
i = 0;
do {
@@ -1359,7 +1393,8 @@ int main(int argc, char *argv[])
n_dif);
if ( (outfile != NULL) && (amb != NULL) ) {
- write_reindexed_stream(infile, outfile, assignments, amb);
+ write_reindexed_stream(infile, outfile, assignments, amb,
+ argc, argv);
} else if ( outfile != NULL ) {
ERROR("Can only write stream with known ambiguity operator.\n");
ERROR("Try again with -w or --operator.\n");
diff --git a/src/geoptimiser.c b/src/geoptimiser.c
index f63ba81b..9b03d1b9 100644
--- a/src/geoptimiser.c
+++ b/src/geoptimiser.c
@@ -85,20 +85,25 @@ static void show_help(const char *s)
" -c, --connected=<rg_coll> Rigid group collection for connected\n"
" ASICs.\n"
" --no-error-maps Do not generate error map PNGs.\n"
+" --stretch-map Generate stretch map PNG (panels distance).\n"
" -x, --min-num-peaks-per-pixel=<num> Minimum number of peaks per pixel.\n"
" Default: 3. \n"
+" -z, --max-num-peaks-per-pixel=<num> Maximum number of peaks per pixel.\n"
+" Default is num_patterns / 10. \n"
" -p, --min-num-pixels-per-conn-group=<num> Minimum number of useful pixels per\n"
" connected group.\n"
-" f Default: 100.\n"
+" Default: 100.\n"
" -l, --most-freq-clen Use only the most frequent camera\n"
" length.\n"
" -s, --individual-dist-offset Use a distance offset for each panel.\n"
" Default: whole-detector offset.\n"
" --no-stretch Do not optimize distance offset.\n"
" Default: distance offset is optimized\n"
+" --no-rotation Do not optimize rotation of connectedd groups.\n"
+" Default: rotation is optimized\n"
" -m --max-peak-dist=<num> Maximum distance between predicted and\n"
" detected peaks (in pixels)\n"
-" Default: half of minimal inter-Bragg distance\n"
+" Default: half of minimal inter-Bragg distance\n"
);
}
@@ -108,11 +113,14 @@ struct geoptimiser_params
char *outfile;
char *geometry_filename;
int min_num_peaks_per_pix;
+ int max_num_peaks_per_pix;
int min_num_pix_per_conn_group;
int only_best_distance;
int nostretch;
+ int norotation;
int individual_coffset;
int error_maps;
+ int stretch_map;
int enforce_cspad_layout;
int no_cspad;
double max_peak_dist;
@@ -808,11 +816,13 @@ static int compute_pixel_displacements(struct image *images, int n_images,
static int compute_avg_pix_displ(struct gpanel *gp, int idx,
- int num_peaks_per_pixel)
+ int num_peaks_per_pixel,
+ int max_peaks_per_pixel)
{
int ret;
- if ( gp->num_pix_displ[idx] >= num_peaks_per_pixel ) {
+ if ( gp->num_pix_displ[idx] >= num_peaks_per_pixel &&
+ gp->num_pix_displ[idx] < max_peaks_per_pixel ) {
ret = fill_avg_pixel_displ(gp, idx);
if ( ret != 0 ) return ret;
@@ -833,12 +843,14 @@ static int compute_avg_pix_displ(struct gpanel *gp, int idx,
static int compute_avg_displacements(struct detector *det,
struct rg_collection *connected,
struct connected_data *conn_data,
- struct gpanel *gpanels)
+ struct gpanel *gpanels,
+ struct geoptimiser_params *gparams)
{
int di, ip, ifs, iss;
- int pix_index, ret;
+ int pix_index, ret, max_peaks;
struct panel *p;
+ max_peaks = 0;
for ( di=0; di<connected->n_rigid_groups; di++ ) {
for ( ip=0; ip<connected->rigid_groups[di]->n_panels; ip++ ) {
@@ -854,9 +866,12 @@ static int compute_avg_displacements(struct detector *det,
for ( ifs=0; ifs<p->w; ifs++ ) {
pix_index = ifs+p->w*iss;
+ if ( max_peaks < gp->num_pix_displ[pix_index] )
+ max_peaks = gp->num_pix_displ[pix_index];
ret = compute_avg_pix_displ(gp, pix_index,
- conn_data[di].num_peaks_per_pixel);
+ conn_data[di].num_peaks_per_pixel,
+ gparams->max_num_peaks_per_pix);
if ( ret != 0 ) return ret;
@@ -865,6 +880,7 @@ static int compute_avg_displacements(struct detector *det,
}
}
+ STATUS("Maximum peaks per pixel is: %d\n", max_peaks);
return 0;
}
@@ -1085,20 +1101,23 @@ static void correct_rotation_and_stretch(struct rg_collection *connected,
/* Calculate corner adjustment
* NB All panels follow the first one */
- if ( ip > 0 && connected->rigid_groups[di]->n_panels == 2 && !gparams->no_cspad ) {
+ if ( ip > 0 ) {
struct panel *p0;
- double delta_x, delta_y, delta;
+ double dx_old, dy_old, dx_new, dy_new;
p0 = connected->rigid_groups[di]->panels[0];
- delta_x = p->cnx - p0->cnx / cs;
- delta_y = p->cny - p0->cny / cs;
+ dx_old = p->cnx - p0->cnx / cs;
+ dy_old = p->cny - p0->cny / cs;
- delta = sqrt(delta_x*delta_x + delta_y*delta_y);
+ dx_new = dx_old*cos(conn_data[di].cang)-
+ dy_old*sin(conn_data[di].cang);
+ dy_new = dx_old*sin(conn_data[di].cang)+
+ dy_old*cos(conn_data[di].cang);
- new_cnx = p0->cnx + delta*p0->fsx;
- new_cny = p0->cny + delta*p0->fsy;
+ new_cnx = p0->cnx + dx_new;
+ new_cny = p0->cny + dy_new;
} else {
@@ -1577,6 +1596,14 @@ static double compute_rotation_and_stretch(struct rg_collection *connected,
free(stretches);
}
+ if ( gparams->norotation ) {
+ STATUS("However, rotation will not be optimized, so the "
+ "rotation angle has been set to 0\n");
+ for ( di=0; di<connected->n_rigid_groups; di++ ) {
+ conn_data[di].cang = 0.0;
+ }
+ }
+
stretch_cf = 1.0;
printf("Computing overall stretch coefficient.\n");
@@ -1642,7 +1669,8 @@ static double compute_rotation_and_stretch(struct rg_collection *connected,
#ifdef CAN_SAVE_TO_PNG
static void draw_panel(struct image *image, cairo_t *cr,
- cairo_matrix_t *basic_m, GdkPixbuf **pixbufs, int i) {
+ cairo_matrix_t *basic_m, GdkPixbuf **pixbufs, int i)
+{
struct panel p = image->det->panels[i];
int w = gdk_pixbuf_get_width(pixbufs[i]);
int h = gdk_pixbuf_get_height(pixbufs[i]);
@@ -1816,8 +1844,119 @@ static int save_data_to_png(char *filename, struct detector *det,
return 0;
}
+
+static int save_stretch_to_png(char *filename, struct detector *det,
+ struct rg_collection *connected,
+ struct connected_data *conn_data)
+
+{
+ struct image im;
+ int i, di, ip;
+ float min_str, max_str;
+ struct rectangle rect;
+ GdkPixbuf *col_scale;
+ cairo_t *cr;
+
+ cairo_status_t r;
+ cairo_surface_t *surf;
+
+ im.det = det;
+ im.bad = NULL;
+ im.dp = malloc(det->n_panels*sizeof(float *));
+ if ( im.dp == NULL ) {
+ ERROR("Failed to allocate data\n");
+ return 1;
+ }
+
+ min_str = 10000;
+ max_str = 0;
+ for ( di=0; di<connected->n_rigid_groups; di++ ) {
+ if (conn_data[di].cstr > max_str) max_str = conn_data[di].cstr;
+ if (conn_data[di].cstr < min_str) min_str = conn_data[di].cstr;
+ }
+
+ i = 0;
+ for ( di=0; di<connected->n_rigid_groups; di++ ) {
+ for ( ip=0; ip<connected->rigid_groups[di]->n_panels; ip++ ) {
+
+ int fs, ss;
+ float val;
+ struct panel *p;
+
+ p = connected->rigid_groups[di]->panels[ip];
+
+ im.dp[i] = calloc(p->w * p->h, sizeof(float));
+ if ( im.dp[i] == NULL ) {
+ ERROR("Failed to allocate data\n");
+ return 1;
+ }
+
+ val = (conn_data[di].cstr - min_str) / (max_str - min_str);
+ val *= 10.0; /* render_panels sets this as max */
+
+ for ( ss=0; ss<p->h; ss++ ) {
+ for ( fs=0; fs<p->w; fs++ ) {
+
+ im.dp[i][fs+p->w*ss] = val;
+
+ }
+ }
+ i++;
+ }
+ }
+
+ get_pixel_extents(im.det, &rect.min_x, &rect.min_y, &rect.max_x,
+ &rect.max_y);
+
+ if (rect.min_x > 0.0) rect.min_x = 0.0;
+ if (rect.max_x < 0.0) rect.max_x = 0.0;
+ if (rect.min_y > 0.0) rect.min_y = 0.0;
+ if (rect.max_y < 0.0) rect.max_y = 0.0;
+
+ rect.width = rect.max_x - rect.min_x;
+ rect.height = rect.max_y - rect.min_y;
+
+ /* Add a thin border */
+ rect.width += 2.0;
+ rect.height += 2.0;
+ surf = cairo_image_surface_create(CAIRO_FORMAT_ARGB32, rect.width + 20,
+ rect.height);
+
+ draw_detector(surf, &im, rect);
+
+ col_scale = render_get_colour_scale(20, rect.height, SCALE_GEOPTIMISER);
+
+ cr = cairo_create(surf);
+ cairo_identity_matrix(cr);
+ cairo_translate(cr, rect.width, 0.0);
+ cairo_rectangle(cr, 0.0, 0.0, 20.0, rect.height);
+ gdk_cairo_set_source_pixbuf(cr, col_scale, 0.0, 0.0);
+ cairo_fill(cr);
+ cairo_destroy(cr);
+
+ for ( i=0; i<det->n_panels; i++ ) {
+ free(im.dp[i]);
+ }
+ free(im.dp);
+
+ r = cairo_surface_write_to_png(surf, filename);
+ if ( r != CAIRO_STATUS_SUCCESS ) return 1;
+
+ return 0;
+}
+
#else /* CAN_SAVE_TO_PNG */
+static int save_stretch_to_png(char *filename, struct detector *det,
+ struct rg_collection *connected,
+ struct connected_data *conn_data)
+{
+ ERROR("WARNING: geoptimiser was compiled without gdk-pixbuf and cairo "
+ "support. Error maps will not be saved.\n");
+ return 0; /* Return "success" so that program continues */
+}
+
+
static int save_data_to_png(char *filename, struct detector *det,
struct gpanel *gpanels)
{
@@ -1902,6 +2041,87 @@ int check_and_enforce_cspad_dist(struct geoptimiser_params *gparams,
}
}
+ return num_errors_found;
+}
+
+
+int check_and_enforce_agipd_dist(struct geoptimiser_params *gparams,
+ struct detector *det)
+{
+ int np = 0;
+ int npp = 0;
+ int num_errors_found = 0;
+
+ double dist_to_check = 66.0;
+ double tol = 0.1;
+
+ for ( np=0; np<det->n_panels; np = np+8 ) {
+
+ double dist2;
+ double cnx2, cny2;
+
+ struct panel *ep = &det->panels[np];
+
+ for ( npp=0; npp<8; npp++ ) {
+
+ struct panel *op = &det->panels[np+npp];
+
+ cnx2 = ep->cnx + npp*dist_to_check*ep->ssx;
+ cny2 = ep->cny + npp*dist_to_check*ep->ssy;
+
+ dist2 = (( cnx2 - op->cnx )*( cnx2 - op->cnx ) +
+ ( cny2 - op->cny )*( cny2 - op->cny ));
+
+ if ( dist2 > (tol*tol)) {
+
+ num_errors_found += 1;
+
+ STATUS("Warning: distance between panels %s and %s "
+ "is outside acceptable margins (Corners are "
+ "more than 0.2 pixels away: %3.2f).\n", ep->name,
+ op->name, sqrt(dist2));
+
+ if ( gparams->enforce_cspad_layout ) {
+
+ double new_op_cx, new_op_cy;
+
+ new_op_cx = ep->cnx + ep->ssx*npp*dist_to_check;
+ new_op_cy = ep->cny + ep->ssy*npp*dist_to_check;
+
+ op->cnx = new_op_cx;
+ op->cny = new_op_cy;
+
+ STATUS("Enforcing distance....\n");
+ }
+
+ }
+
+ if ( ep->fsx != op->fsx || ep->ssx != op->ssx ||
+ ep->fsy != op->fsy || ep->ssx != op->ssx ) {
+
+ num_errors_found += 1;
+
+ STATUS("Warning: relative orientation between panels "
+ "%s and %s is incorrect.\n", ep->name, op->name);
+
+ if ( gparams->enforce_cspad_layout ) {
+
+ STATUS("Enforcing relative orientation....\n");
+
+ op->fsx = ep->fsx;
+ op->ssx = ep->ssx;
+ op->fsy = ep->fsy;
+ op->ssy = ep->ssy;
+
+ op->xfs = ep->xfs;
+ op->xss = ep->xss;
+ op->yfs = ep->yfs;
+ op->yss = ep->yss;
+ }
+
+ }
+ }
+ }
return num_errors_found;
}
@@ -2076,9 +2296,9 @@ int optimize_geometry(struct geoptimiser_params *gparams,
double avg_res;
double clen_to_use;
- // for angles and stretch calculation use
- // only pixels which are distco*size_panel
- // away
+ /* for angles and stretch calculation use
+ * only pixels which are distco*size_panel
+ * away */
double dist_coeff_for_rot_str = 0.2;
double total_error;
@@ -2099,6 +2319,7 @@ int optimize_geometry(struct geoptimiser_params *gparams,
"accurate estimation of position/orientation: %i\n",
gparams->min_num_pix_per_conn_group);
+ /* CS-PAD */
if ( (det->n_panels == 64) && !gparams->no_cspad ) {
int num_errors = 0;
@@ -2142,6 +2363,51 @@ int optimize_geometry(struct geoptimiser_params *gparams,
}
}
+ /* AGIPD */
+ if ( (det->n_panels == 128) && det->panels[0].res > 4999 &&
+ det->panels[0].res < 5001 && !gparams->no_cspad ) {
+
+ int num_errors = 0;
+
+ STATUS("It looks like the detector is an AGIPD. "
+ "Checking relative distance and orientation of "
+ "connected ASICS.\n");
+ STATUS("If the detector is not an AGIPD , please rerun the "
+ "program with the --no-agipd option.\n");
+
+ STATUS("Enforcing AGIPD layout...\n");
+ num_errors = check_and_enforce_agipd_dist(gparams, det);
+
+ if ( gparams->enforce_cspad_layout ) {
+
+ int geom_wr;
+
+ STATUS("Saving geometry with enforced AGIPD layout.\n"
+ "Please restart geometry optimization using the "
+ "optimized geometry from this run as input "
+ "geometry file.\n");
+ geom_wr = write_detector_geometry_2(
+ gparams->geometry_filename,
+ gparams->outfile, det,
+ gparams->command_line, 1);
+ if ( geom_wr != 0 ) {
+ ERROR("Error in writing output geometry file.\n");
+ return 1;
+ }
+ STATUS("All done!\n");
+ return 0;
+ }
+
+ if ( !gparams->enforce_cspad_layout && num_errors > 0 ) {
+ ERROR("Relative distance and orientation of connected "
+ "ASICS do not respect the AGIPD layout.\n"
+ "Geometry optimization cannot continue.\n"
+ "Please rerun the program with the "
+ "--enforce-agipd-layout option.\n");
+ return 1;
+ }
+ }
+
images = read_patterns_from_stream(gparams->infile, det, &n_images);
if ( (n_images < 1) || (images == NULL) ) {
ERROR("Error reading stream file\n");
@@ -2164,6 +2430,14 @@ int optimize_geometry(struct geoptimiser_params *gparams,
avg_cell);
if ( clen_to_use < 0.0 ) return 1;
+ if ( gparams->max_num_peaks_per_pix == 0 && n_images > 100 ) {
+ gparams->max_num_peaks_per_pix = n_images / 10;
+ } else if ( gparams->max_num_peaks_per_pix == 0 ) {
+ gparams->max_num_peaks_per_pix = n_images;
+ }
+ STATUS("Maximum number of measurements for a pixel to be included in "
+ "the refinement: %i\n", gparams->max_num_peaks_per_pix);
+
gpanels = calloc(det->n_panels, sizeof(struct gpanel));
if ( gpanels == NULL ) {
ERROR("Failed to allocate panels.\n");
@@ -2206,7 +2480,7 @@ int optimize_geometry(struct geoptimiser_params *gparams,
adjust_min_peaks_per_conn(connected, gpanels, det, gparams, conn_data);
- if ( compute_avg_displacements(det, connected, conn_data, gpanels) ) {
+ if ( compute_avg_displacements(det, connected, conn_data, gpanels, gparams) ) {
free(conn_data);
free(images);
return 1;
@@ -2235,25 +2509,53 @@ int optimize_geometry(struct geoptimiser_params *gparams,
}
- stretch_coeff = compute_rotation_and_stretch(connected, conn_data,
- det, gpanels,
+ if ( !gparams->nostretch || !gparams->norotation ) {
+
+ stretch_coeff = compute_rotation_and_stretch(connected,
+ conn_data, det, gpanels,
dist_coeff_for_rot_str,
gparams);
- if ( stretch_coeff < 0.0 ) {
- free(conn_data);
- return 1;
- }
+ if ( stretch_coeff < 0.0 ) {
+ free(conn_data);
+ return 1;
+ }
- ret = compute_rot_stretch_for_empty_panels(quadrants, connected,
- gparams->min_num_pix_per_conn_group, conn_data);
- if ( ret ) {
- free(conn_data);
- return 1;
- }
+ ret = compute_rot_stretch_for_empty_panels(quadrants, connected,
+ gparams->min_num_pix_per_conn_group,
+ conn_data);
+ if ( ret ) {
+ free(conn_data);
+ return 1;
+ }
+
+ if ( gparams->stretch_map ) {
+
+#ifdef HAVE_SAVE_TO_PNG
+
+ STATUS("Saving stretch map - for out-of-plane rotations.\n");
- correct_rotation_and_stretch(connected, det, gpanels, conn_data,
- clen_to_use, stretch_coeff,
- gparams);
+ ret = save_stretch_to_png("stretch_co.png", det,
+ connected, conn_data);
+
+ if ( ret ) {
+ ERROR("Error while writing data to file.\n");
+ free(conn_data);
+ return 1;
+ }
+
+#else /* HAVE_SAVE_TO_PNG */
+
+ STATUS("ERROR: geoptimiser was compiled without GTK and "
+ "cairo support.\n Stretch map will not be saved.\n");
+
+#endif /* HAVE_SAVE_TO_PNG */
+
+ }
+
+ correct_rotation_and_stretch(connected, det, gpanels,
+ conn_data, clen_to_use,
+ stretch_coeff, gparams);
+ }
ret = compute_shift(connected, conn_data, det, gparams, gpanels);
if ( ret != 0 ) {
@@ -2341,13 +2643,16 @@ int main(int argc, char *argv[])
gparams->infile = NULL;
gparams->geometry_filename = NULL;
gparams->min_num_peaks_per_pix = 3;
+ gparams->max_num_peaks_per_pix = 0;
gparams->min_num_pix_per_conn_group = 100;
gparams->only_best_distance = 0;
gparams->enforce_cspad_layout = 0;
gparams->nostretch = 0;
+ gparams->norotation = 0;
gparams->individual_coffset = 0;
gparams->no_cspad = 0;
gparams->error_maps = 1;
+ gparams->stretch_map = 0;
gparams->max_peak_dist = 0.0;
const struct option longopts[] = {
@@ -2361,6 +2666,7 @@ int main(int argc, char *argv[])
{"quadrants", 1, NULL, 'q'},
{"connected", 1, NULL, 'c'},
{"min-num-peaks-per-pixel", 1, NULL, 'x'},
+ {"max-num-peaks-per-pixel", 0, NULL, 'z'},
{"min-num-pixels-per-conn-group", 1, NULL, 'p'},
{"most-few-clen", 0, NULL, 'l'},
{"max-peak-dist", 1, NULL, 'm'},
@@ -2368,9 +2674,13 @@ int main(int argc, char *argv[])
/* Long-only options with no arguments */
{"no-stretch", 0, &gparams->nostretch, 1},
+ {"no-rotation", 0, &gparams->norotation, 1},
{"no-error-maps", 0, &gparams->error_maps, 0},
+ {"stretch-map", 0, &gparams->stretch_map, 1},
{"enforce-cspad-layout", 0, &gparams->enforce_cspad_layout, 1},
+ {"enforce-agipd-layout", 0, &gparams->enforce_cspad_layout, 1},
{"no-cspad", 0, &gparams->no_cspad, 1},
+ {"no-agipd", 0, &gparams->no_cspad, 1},
/* Long-only options with arguments */
{"min-num-peaks-per-panel", 1, NULL, 11},
@@ -2380,7 +2690,7 @@ int main(int argc, char *argv[])
};
/* Short options */
- while ((c = getopt_long(argc, argv, "ho:i:g:q:c:o:x:p:lsm:",
+ while ((c = getopt_long(argc, argv, "ho:i:g:q:c:o:x:z:p:lsm:",
longopts, NULL)) != -1) {
switch (c) {
@@ -2425,6 +2735,10 @@ int main(int argc, char *argv[])
gparams->min_num_peaks_per_pix = atoi(optarg);
break;
+ case 'z' :
+ gparams->max_num_peaks_per_pix = atoi(optarg);
+ break;
+
case 'p' :
gparams->min_num_pix_per_conn_group = atoi(optarg);
break;
diff --git a/src/im-sandbox.c b/src/im-sandbox.c
index e01fcb74..74467529 100644
--- a/src/im-sandbox.c
+++ b/src/im-sandbox.c
@@ -196,7 +196,7 @@ static struct filename_plus_event *get_pattern(FILE *fh, int config_basename,
line = tmp;
}
- scan_check = sscanf(line, "%s %s", filename_buf, event_buf );
+ scan_check = sscanf(line, "%s %s", filename_buf, event_buf);
len = strlen(prefix)+strlen(filename_buf)+1;
@@ -261,9 +261,15 @@ static struct filename_plus_event *get_pattern(FILE *fh, int config_basename,
struct event *ev_to_add;
ev_to_add = get_event_from_event_string(event_buf);
- append_event_to_event_list(ev_list, ev_to_add);
- free_event(ev_to_add);
- event_index = 0;
+ if ( ev_to_add == NULL ) {
+ ERROR("Bad event descriptor: '%s'\n",
+ event_buf);
+ } else {
+ append_event_to_event_list(ev_list,
+ ev_to_add);
+ free_event(ev_to_add);
+ event_index = 0;
+ }
}
diff --git a/src/indexamajig.c b/src/indexamajig.c
index bac1ee04..38d83d58 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -123,11 +123,11 @@ static void show_help(const char *s)
" --no-check-cell Don't check lattice parameters against input cell\n"
" --no-cell-combinations\n"
" Don't use axis combinations when checking cell\n"
-" --no-multi Don't repeat indexing to index multiple hits\n"
+" --multi Repeat indexing to index multiple hits\n"
" --no-retry Don't repeat indexing to increase indexing rate\n"
" --no-refine Skip the prediction refinement step\n"
-" --check-peaks Check that most of the peaks can be accounted for\n"
-" by the indexing solution\n"
+" --no-check-peaks Don't check that most of the peaks can be accounted\n"
+" for by the indexing solution\n"
"\n"
" --taketwo-member-threshold\n"
" Minimum number of members in network\n"
@@ -238,8 +238,8 @@ int main(int argc, char *argv[])
int if_refine = 1;
int if_nocomb = 0;
int if_nocheck = 0;
- int if_peaks = 0;
- int if_multi = 1;
+ int if_peaks = 1;
+ int if_multi = 0;
int if_retry = 1;
int serial_start = 1;
@@ -342,8 +342,11 @@ int main(int argc, char *argv[])
{"no-check-cell", 0, &if_nocheck, 1},
{"no-cell-check", 0, &if_nocheck, 1},
{"check-peaks", 0, &if_peaks, 1},
+ {"no-check-peaks", 0, &if_peaks, 0},
{"no-retry", 0, &if_retry, 0},
+ {"retry", 0, &if_retry, 1},
{"no-multi", 0, &if_multi, 0},
+ {"multi", 0, &if_multi, 1},
{"overpredict", 0, &iargs.overpredict, 1},
/* Long-only options which don't actually do anything */
@@ -851,6 +854,8 @@ int main(int argc, char *argv[])
if ( have_push_res && !(iargs.int_meth & INTEGRATION_RESCUT) ) {
ERROR("WARNING: You used --push-res, but not -rescut, "
"therefore --push-res will have no effect.\n");
+ ERROR("WARNING: Add --integration=rings-rescut or "
+ "--integration=prof2d-rescut.\n");
}
/* Parse unit cell tolerance */
@@ -955,15 +960,6 @@ int main(int argc, char *argv[])
}
- /* Open output stream */
- st = open_stream_for_write_3(outfile, geom_filename, iargs.cell,
- argc, argv);
- if ( st == NULL ) {
- ERROR("Failed to open stream '%s'\n", outfile);
- return 1;
- }
- free(outfile);
-
if ( indm_str == NULL ) {
STATUS("No indexing methods specified. I will try to ");
@@ -1035,10 +1031,19 @@ int main(int argc, char *argv[])
ERROR("Failed to set up indexing system\n");
return 1;
}
- free(indm_str);
}
+ /* Open output stream */
+ st = open_stream_for_write_4(outfile, geom_filename, iargs.cell,
+ argc, argv, indm_str);
+ if ( st == NULL ) {
+ ERROR("Failed to open stream '%s'\n", outfile);
+ return 1;
+ }
+ free(outfile);
+ free(indm_str);
+
gsl_set_error_handler_off();
create_sandbox(&iargs, n_proc, prefix, config_basename, fh,
diff --git a/src/merge.c b/src/merge.c
index 17e57c90..14460953 100644
--- a/src/merge.c
+++ b/src/merge.c
@@ -50,7 +50,7 @@
/* Minimum partiality of a reflection for it to be merged */
-#define MIN_PART_MERGE (0.05)
+#define MIN_PART_MERGE (0.3)
struct merge_queue_args
@@ -59,7 +59,6 @@ struct merge_queue_args
pthread_rwlock_t full_lock;
Crystal **crystals;
int n_started;
- PartialityModel pmodel;
double push_res;
int use_weak;
long long int n_reflections;
@@ -159,7 +158,6 @@ static void run_merge_job(void *vwargs, int cookie)
signed int h, k, l;
double mean, sumweight, M2, temp, delta, R;
double corr, res, w;
- //double esd;
if ( get_partiality(refl) < MIN_PART_MERGE ) continue;
@@ -189,7 +187,7 @@ static void run_merge_job(void *vwargs, int cookie)
}
/* Total (multiplicative) correction factor */
- corr = exp(-G) * exp(B*res*res) * get_lorentz(refl)
+ corr = 1.0/G * exp(B*res*res) * get_lorentz(refl)
/ get_partiality(refl);
if ( isnan(corr) ) {
ERROR("NaN corr:\n");
@@ -200,7 +198,7 @@ static void run_merge_job(void *vwargs, int cookie)
continue;
}
- //esd = get_esd_intensity(refl) * corr;
+ /* Reflections count less the more they have to be scaled up */
w = 1.0;
/* Running mean and variance calculation */
@@ -229,7 +227,7 @@ static void finalise_merge_job(void *vqargs, void *vwargs)
RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
- PartialityModel pmodel, int min_meas,
+ int min_meas,
double push_res, int use_weak)
{
RefList *full;
@@ -245,7 +243,6 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
qargs.full = full;
qargs.n_started = 0;
qargs.crystals = crystals;
- qargs.pmodel = pmodel;
qargs.push_res = push_res;
qargs.use_weak = use_weak;
qargs.n_reflections = 0;
@@ -285,3 +282,120 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
reflist_free(full);
return full2;
}
+
+
+double residual(Crystal *cr, const RefList *full, int free,
+ int *pn_used, const char *filename)
+{
+ Reflection *refl;
+ RefListIterator *iter;
+ int n_used = 0;
+ double num = 0.0;
+ double den = 0.0;
+ double G = crystal_get_osf(cr);
+ double B = crystal_get_Bfac(cr);
+ UnitCell *cell = crystal_get_cell(cr);
+
+ for ( refl = first_refl(crystal_get_reflections(cr), &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ double p, w, corr, res;
+ signed int h, k, l;
+ Reflection *match;
+ double I_full;
+ double int1, int2;
+
+ if ( free && !get_flag(refl) ) continue;
+
+ get_indices(refl, &h, &k, &l);
+ res = resolution(cell, h, k, l);
+ match = find_refl(full, h, k, l);
+ if ( match == NULL ) continue;
+ I_full = get_intensity(match);
+
+ if ( get_redundancy(match) < 2 ) continue;
+
+ p = get_partiality(refl);
+ //if ( p < 0.2 ) continue;
+
+ corr = get_lorentz(refl) / (G * exp(-B*res*res));
+ int1 = get_intensity(refl) * corr;
+ int2 = p*I_full;
+ w = p;
+
+ num += fabs(int1 - int2) * w;
+ den += fabs(int1 + int2) * w/2.0;
+
+ n_used++;
+
+ }
+
+ if ( pn_used != NULL ) *pn_used = n_used;
+ return num/(den*sqrt(2));
+}
+
+
+double log_residual(Crystal *cr, const RefList *full, int free,
+ int *pn_used, const char *filename)
+{
+ double dev = 0.0;
+ double G, B;
+ Reflection *refl;
+ RefListIterator *iter;
+ int n_used = 0;
+ FILE *fh = NULL;
+
+ G = crystal_get_osf(cr);
+ B = crystal_get_Bfac(cr);
+ if ( filename != NULL ) {
+ fh = fopen(filename, "a");
+ if ( fh == NULL ) {
+ ERROR("Failed to open '%s'\n", filename);
+ }
+ }
+
+ for ( refl = first_refl(crystal_get_reflections(cr), &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ double p, L, s, w;
+ signed int h, k, l;
+ Reflection *match;
+ double esd, I_full, I_partial;
+ double fx;
+
+ if ( free && !get_flag(refl) ) continue;
+
+ get_indices(refl, &h, &k, &l);
+ match = find_refl(full, h, k, l);
+ if ( match == NULL ) continue;
+
+ p = get_partiality(refl);
+ L = get_lorentz(refl);
+ I_partial = get_intensity(refl);
+ I_full = get_intensity(match);
+ esd = get_esd_intensity(refl);
+ s = resolution(crystal_get_cell(cr), h, k, l);
+
+ if ( I_partial <= 3.0*esd ) continue; /* Also because of log */
+ if ( get_redundancy(match) < 2 ) continue;
+ if ( I_full <= 0 ) continue; /* Because log */
+ if ( p <= 0.0 ) continue; /* Because of log */
+
+ fx = log(G) - B*s*s + log(p) + log(I_full) - log(I_partial) - log(L);
+ w = 1.0;
+ dev += w*fx*fx;
+
+ if ( fh != NULL ) {
+ fprintf(fh, "%4i %4i %4i %e %e\n",
+ h, k, l, s, dev);
+ }
+
+ }
+
+ if ( fh != NULL ) fclose(fh);
+
+ if ( pn_used != NULL ) *pn_used = n_used;
+ return dev;
+}
diff --git a/src/merge.h b/src/merge.h
index 1911584f..b33afdea 100644
--- a/src/merge.h
+++ b/src/merge.h
@@ -40,7 +40,12 @@
#include "geometry.h"
extern RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
- PartialityModel pmodel, int min_meas,
- double push_res, int use_weak);
+ int min_meas, double push_res, int use_weak);
+
+extern double residual(Crystal *cr, const RefList *full, int free,
+ int *pn_used, const char *filename);
+
+extern double log_residual(Crystal *cr, const RefList *full, int free,
+ int *pn_used, const char *filename);
#endif /* MERGE */
diff --git a/src/partial_sim.c b/src/partial_sim.c
index 5926043f..218780fd 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -385,7 +385,7 @@ static void run_job(void *vwargs, int cookie)
reflections = predict_to_res(cr, largest_q(&wargs->image));
crystal_set_reflections(cr, reflections);
- calculate_partialities(cr, PMODEL_SCSPHERE);
+ calculate_partialities(cr, PMODEL_XSPHERE);
for ( i=0; i<NBINS; i++ ) {
wargs->n_ref[i] = 0;
@@ -861,7 +861,7 @@ int main(int argc, char *argv[])
STATUS(" Background: none (no image "
"output)\n");
}
- STATUS(" Partiality model: scsphere (hardcoded)\n");
+ STATUS(" Partiality model: xsphere (hardcoded)\n");
STATUS(" Noise standard deviation: %.2f detector units\n",
noise_stddev);
if ( random_intensities ) {
@@ -873,6 +873,7 @@ int main(int argc, char *argv[])
input_file, sym_str);
}
STATUS(" Max error in cell components: %.2f %%\n", cnoise);
+ STATUS("Scale factor standard deviation: %.2f\n", osf_stddev);
if ( random_intensities ) {
full = reflist_new();
diff --git a/src/partialator.c b/src/partialator.c
index 2540521d..e152db7c 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -3,11 +3,11 @@
*
* Scaling and post refinement for coherent nanocrystallography
*
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2015 Thomas White <taw@physics.org>
+ * 2010-2018 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -40,6 +40,7 @@
#include <assert.h>
#include <pthread.h>
#include <gsl/gsl_errno.h>
+#include <sys/stat.h>
#include <image.h>
#include <utils.h>
@@ -189,7 +190,7 @@ static void write_split(Crystal **crystals, int n_crystals, const char *outfile,
}
snprintf(tmp, 1024, "%s1", outfile);
split = merge_intensities(crystals1, n_crystals1, nthreads,
- pmodel, min_measurements, push_res, 1);
+ min_measurements, push_res, 1);
if ( split == NULL ) {
ERROR("Not enough crystals for two way split!\n");
@@ -201,7 +202,7 @@ static void write_split(Crystal **crystals, int n_crystals, const char *outfile,
reflist_free(split);
snprintf(tmp, 1024, "%s2", outfile);
split = merge_intensities(crystals2, n_crystals2, nthreads,
- pmodel, min_measurements, push_res, 1);
+ min_measurements, push_res, 1);
STATUS("and %s\n", tmp);
write_reflist_2(tmp, split, sym);
reflist_free(split);
@@ -288,7 +289,7 @@ static void write_custom_split(struct custom_split *csplit, int dsn,
STATUS("Writing dataset '%s' to %s (%i crystals)\n",
csplit->dataset_names[dsn], tmp, n_crystalsn);
split = merge_intensities(crystalsn, n_crystalsn, nthreads,
- pmodel, min_measurements, push_res, 1);
+ min_measurements, push_res, 1);
write_reflist_2(tmp, split, sym);
reflist_free(split);
@@ -315,6 +316,7 @@ static void show_help(const char *s)
" --stop-after=<n> Stop after merging <n> crystals.\n"
" -n, --iterations=<n> Run <n> cycles of scaling and post-refinement.\n"
" --no-scale Disable scale factor (G, B) refinement.\n"
+" --no-Bscale Disable B factor scaling.\n"
" --no-pr Disable orientation/physics refinement.\n"
" -m, --model=<model> Specify partiality model.\n"
" --min-measurements=<n> Minimum number of measurements to require.\n"
@@ -324,7 +326,12 @@ static void show_help(const char *s)
" -j <n> Run <n> analyses in parallel.\n"
" --no-free Disable cross-validation (testing only).\n"
" --custom-split List of files for custom dataset splitting.\n"
-" --max-rel-B Maximum allowable relative |B| factor.\n");
+" --max-rel-B Maximum allowable relative |B| factor.\n"
+" --no-logs Do not write extensive log files.\n"
+" -w <pg> Apparent point group for resolving ambiguities.\n"
+" --operator=<op> Indexing ambiguity operator for resolving.\n"
+" --force-bandwidth=<n> Set all bandwidths to <n> (fraction).\n"
+" --force-radius=<n> Set all profile radii to <n> nm^-1.\n");
}
@@ -578,7 +585,7 @@ static int set_initial_params(Crystal *cr, FILE *fh)
} else {
- crystal_set_osf(cr, 0.0);
+ crystal_set_osf(cr, 1.0);
crystal_set_Bfac(cr, 0.0);
}
@@ -603,13 +610,22 @@ static void select_free_reflections(RefList *list, gsl_rng *rng)
static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
- int fr)
+ int fr, signed int inum)
{
Reflection *refl;
RefListIterator *iter;
double G = crystal_get_osf(cr);
double B = crystal_get_Bfac(cr);
UnitCell *cell = crystal_get_cell(cr);
+ char ins[5];
+
+ if ( inum >= 0 ) {
+ snprintf(ins, 4, "%i", inum);
+ } else {
+ ins[0] = 'F';
+ ins[1] = '\0';
+ }
+
for ( refl = first_refl(list, &iter);
refl != NULL;
@@ -622,6 +638,9 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
if ( !get_flag(refl) ) continue; /* Not free-flagged */
+ /* Strong reflections only */
+ if ( get_intensity(refl) < 3.0*get_esd_intensity(refl) ) continue;
+
get_indices(refl, &h, &k, &l);
res = resolution(cell, h, k, l);
if ( 2.0*res > crystal_get_resolution_limit(cr) ) continue;
@@ -629,46 +648,55 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
match = find_refl(full, h, k, l);
if ( match == NULL ) continue;
+ /* Don't calculate pobs if reference reflection is weak */
+ if ( fabs(get_intensity(match)) / get_esd_intensity(match) < 3.0 ) continue;
+
/* Calculated partiality */
pcalc = get_partiality(refl);
/* Observed partiality */
- corr = exp(-G) * exp(B*res*res) * get_lorentz(refl);
+ corr = G * exp(B*res*res) * get_lorentz(refl);
Ipart = get_intensity(refl) * corr;
pobs = Ipart / get_intensity(match);
- fprintf(fh, "%5i %4i %4i %4i %8.4f %8.3f %8.3f\n",
- fr, h, k, l, 2*res/1e9, pcalc, pobs);
+ fprintf(fh, "%5i %4i %4i %4i %e %e %8.3f %8.3f %s\n",
+ fr, h, k, l, 2*res, Ipart, pcalc, pobs, ins);
}
}
static void write_pgraph(RefList *full, Crystal **crystals, int n_crystals,
- int iter)
+ signed int iter, const char *suff)
{
FILE *fh;
char tmp[256];
int i;
- snprintf(tmp, 256, "pgraph-iter%i.dat", iter);
+ snprintf(tmp, 256, "pr-logs/pgraph%s.dat", suff);
+
+ if ( iter == 0 ) {
+ fh = fopen(tmp, "w");
+ } else {
+ fh = fopen(tmp, "a");
+ }
- fh = fopen(tmp, "w");
if ( fh == NULL ) {
ERROR("Failed to open '%s'\n", tmp);
return;
}
- fprintf(fh, " fr h k l 1/d(nm) pcalc pobs\n");
+ if ( iter == 0 ) {
+ fprintf(fh, " Crystal h k l 1/d(m) Ipart pcalc pobs iteration\n");
+ }
for ( i=0; i<n_crystals; i++ ) {
if ( crystal_get_user_flag(crystals[i]) != 0 ) continue;
write_to_pgraph(fh, crystal_get_reflections(crystals[i]), full,
- crystals[i], i);
+ crystals[i], i, iter);
}
fclose(fh);
-
}
@@ -689,6 +717,7 @@ static void all_residuals(Crystal **crystals, int n_crystals, RefList *full,
if ( crystal_get_user_flag(crystals[i]) ) continue;
+ /* Scaling should have been done right before calling this */
r = residual(crystals[i], full, 0, NULL, NULL);
free_r = residual(crystals[i], full, 1, NULL, NULL);
log_r = log_residual(crystals[i], full, 0, NULL, NULL);
@@ -719,32 +748,82 @@ static void show_all_residuals(Crystal **crystals, int n_crystals,
}
-static void dump_parameters(const char *filename, Crystal **crystals,
- int n_crystals)
+struct log_qargs
{
- FILE *fh;
- fh = fopen(filename, "w");
- if ( fh == NULL ) {
- ERROR("Couldn't open partialator.params!\n");
- } else {
- fprintf(fh, " cr OSF relB div"
- " flag filename event\n");
- int i;
- for ( i=0; i<n_crystals; i++ ) {
- struct image *img;
- char *evt_str;
- img = crystal_get_image(crystals[i]);
- evt_str = get_event_string(img->event);
- fprintf(fh, "%4i %10.5f %10.2f %8.5e %-25s %s %s\n",
- i, crystal_get_osf(crystals[i]),
- crystal_get_Bfac(crystals[i])*1e20,
- crystal_get_image(crystals[i])->div,
- str_prflag(crystal_get_user_flag(crystals[i])),
- img->filename, evt_str);
- free(evt_str);
- }
- fclose(fh);
- }
+ int iter;
+ int next;
+ Crystal **crystals;
+ int n_crystals;
+ RefList *full;
+ int scaleflags;
+};
+
+
+struct log_args
+{
+ Crystal *cr;
+ RefList *full;
+ int scaleflags;
+ int iter;
+ int cnum;
+};
+
+
+static void *get_log_task(void *vp)
+{
+ struct log_qargs *qargs = vp;
+ struct log_args *task;
+
+ if ( qargs->next >= qargs->n_crystals ) return NULL;
+
+ task = malloc(sizeof(struct log_args));
+ if ( task == NULL ) return NULL;
+
+ task->cr = qargs->crystals[qargs->next];
+ task->full = qargs->full;
+ task->iter = qargs->iter;
+ task->cnum = qargs->next;
+ task->scaleflags = qargs->scaleflags;
+
+ qargs->next += 20;
+ return task;
+}
+
+
+static void write_logs(void *vp, int cookie)
+{
+ struct log_args *args = vp;
+ write_specgraph(args->cr, args->full, args->iter, args->cnum);
+ write_gridscan(args->cr, args->full, args->iter, args->cnum,
+ args->scaleflags);
+ write_test_logs(args->cr, args->full, args->iter, args->cnum);
+}
+
+
+static void done_log(void *qargs, void *vp)
+{
+ struct log_args *task = vp;
+ free(task);
+}
+
+
+static void write_logs_parallel(Crystal **crystals, int n_crystals,
+ RefList *full, int iter, int n_threads,
+ int scaleflags)
+{
+ struct log_qargs qargs;
+
+ qargs.iter = iter;
+ qargs.next = 0;
+ qargs.full = full;
+ qargs.crystals = crystals;
+ qargs.n_crystals = n_crystals;
+ qargs.scaleflags = scaleflags;
+
+ STATUS("Writing logs...\n");
+ run_threads(n_threads, write_logs, get_log_task, done_log, &qargs,
+ n_crystals/20, 0, 0, 0);
+ STATUS("Done.\n");
}
@@ -755,9 +834,10 @@ int main(int argc, char *argv[])
char *outfile = NULL;
char *sym_str = NULL;
SymOpList *sym;
+ SymOpList *amb;
+ SymOpList *w_sym;
int nthreads = 1;
int i;
- struct image *images;
int n_iter = 10;
RefList *full;
int n_images = 0;
@@ -765,11 +845,12 @@ int main(int argc, char *argv[])
int n_crystals_seen = 0;
char cmdline[1024];
int no_scale = 0;
+ int no_Bscale = 0;
int no_pr = 0;
Stream *st;
Crystal **crystals;
char *pmodel_str = NULL;
- PartialityModel pmodel = PMODEL_SCSPHERE;
+ PartialityModel pmodel = PMODEL_XSPHERE;
int min_measurements = 2;
char *rval;
int polarisation = 1;
@@ -785,6 +866,15 @@ int main(int argc, char *argv[])
char *csplit_fn = NULL;
struct custom_split *csplit = NULL;
double max_B = 1e-18;
+ char *rfile = NULL;
+ RefList *reference = NULL;
+ int no_logs = 0;
+ char *w_sym_str = NULL;
+ char *operator = NULL;
+ double force_bandwidth = -1.0;
+ double force_radius = -1.0;
+ char *audit_info;
+ int scaleflags = 0;
/* Long options */
const struct option longopts[] = {
@@ -805,10 +895,15 @@ int main(int argc, char *argv[])
{"push-res", 1, NULL, 5},
{"res-push", 1, NULL, 5}, /* compat */
{"custom-split", 1, NULL, 6},
- {"max-rel-B" , 1, NULL, 7},
- {"max-rel-b" , 1, NULL, 7}, /* compat */
+ {"max-rel-B", 1, NULL, 7},
+ {"max-rel-b", 1, NULL, 7}, /* compat */
+ {"reference", 1, NULL, 8}, /* ssshhh! */
+ {"operator", 1, NULL, 9},
+ {"force-bandwidth", 1, NULL, 10},
+ {"force-radius", 1, NULL, 11},
{"no-scale", 0, &no_scale, 1},
+ {"no-Bscale", 0, &no_Bscale, 1},
{"no-pr", 0, &no_pr, 1},
{"no-polarisation", 0, &polarisation, 0},
{"no-polarization", 0, &polarisation, 0}, /* compat */
@@ -816,6 +911,7 @@ int main(int argc, char *argv[])
{"polarization", 0, &polarisation, 1}, /* compat */
{"no-free", 0, &no_free, 1},
{"output-every-cycle", 0, &output_everycycle, 1},
+ {"no-logs", 0, &no_logs, 1},
{0, 0, NULL, 0}
};
@@ -827,7 +923,7 @@ int main(int argc, char *argv[])
}
/* Short options */
- while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:j:m:",
+ while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:j:m:w:",
longopts, NULL)) != -1)
{
@@ -886,6 +982,10 @@ int main(int argc, char *argv[])
pmodel_str = strdup(optarg);
break;
+ case 'w' :
+ w_sym_str = strdup(optarg);
+ break;
+
case 2 :
errno = 0;
min_measurements = strtod(optarg, &rval);
@@ -932,6 +1032,33 @@ int main(int argc, char *argv[])
max_B = max_B * 1e-20;
break;
+ case 8 :
+ rfile = strdup(optarg);
+ break;
+
+ case 9 :
+ operator = strdup(optarg);
+ break;
+
+ case 10 :
+ errno = 0;
+ force_bandwidth = strtod(optarg, &rval);
+ if ( *rval != '\0' ) {
+ ERROR("Invalid value for --force-bandwidth.\n");
+ return 1;
+ }
+ break;
+
+ case 11 :
+ errno = 0;
+ force_radius = strtod(optarg, &rval);
+ if ( *rval != '\0' ) {
+ ERROR("Invalid value for --force-radius.\n");
+ return 1;
+ }
+ force_radius *= 1e9;
+ break;
+
case 0 :
break;
@@ -971,13 +1098,49 @@ int main(int argc, char *argv[])
free(sym_str);
if ( sym == NULL ) return 1;
+ if ( (w_sym_str != NULL) && (operator != NULL) ) {
+ ERROR("Specify the apparent symmetry (-w) or the operator, "
+ "not both.\n");
+ return 1;
+ }
+
+ if ( w_sym_str == NULL ) {
+ w_sym = NULL;
+ amb = NULL;
+ } else {
+ pointgroup_warning(w_sym_str);
+ w_sym = get_pointgroup(w_sym_str);
+ free(w_sym_str);
+ if ( w_sym == NULL ) return 1;
+ amb = get_ambiguities(w_sym, sym);
+ if ( amb == NULL ) {
+ ERROR("Couldn't find ambiguity operator.\n");
+ ERROR("Check that your values for -y and -w are "
+ "correct.\n");
+ return 1;
+ }
+
+ }
+
+ if ( operator ) {
+ amb = parse_symmetry_operations(operator);
+ if ( amb == NULL ) return 1;
+ set_symmetry_name(amb, "Ambiguity");
+ }
+
+ if ( amb != NULL ) {
+ STATUS("Indexing ambiguity resolution enabled. "
+ "The ambiguity operation(s) are:\n");
+ describe_symmetry(amb);
+ /* In contrast to ambigator, partialator can deal with multiple
+ * ambiguities at once */
+ }
+
if ( pmodel_str != NULL ) {
if ( strcmp(pmodel_str, "unity") == 0 ) {
pmodel = PMODEL_UNITY;
- } else if ( strcmp(pmodel_str, "scgaussian") == 0 ) {
- pmodel = PMODEL_SCGAUSSIAN;
- } else if ( strcmp(pmodel_str, "scsphere") == 0 ) {
- pmodel = PMODEL_SCSPHERE;
+ } else if ( strcmp(pmodel_str, "xsphere") == 0 ) {
+ pmodel = PMODEL_XSPHERE;
} else if ( strcmp(pmodel_str, "random") == 0 ) {
pmodel = PMODEL_RANDOM;
} else {
@@ -988,7 +1151,25 @@ int main(int argc, char *argv[])
if ( (pmodel == PMODEL_UNITY) && !no_pr ) {
no_pr = 1;
- STATUS("Setting --no-pr because we are not modelling partialities (--model=unity).\n");
+ STATUS("Setting --no-pr because we are not modelling "
+ "partialities (--model=unity).\n");
+ }
+
+ if ( no_Bscale ) {
+ scaleflags |= SCALE_NO_B;
+ }
+
+ if ( !no_logs ) {
+ int r = mkdir("pr-logs", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+ if ( r ) {
+ if ( errno == EEXIST ) {
+ ERROR("WARNING: pr-logs folder already exists. "
+ "Beware of mixing old and new log files!\n");
+ } else {
+ ERROR("Failed to create pr-logs folder.\n");
+ return 1;
+ }
+ }
}
/* Read the custom split list (if applicable) */
@@ -1001,6 +1182,20 @@ int main(int argc, char *argv[])
free(csplit_fn);
}
+ if ( rfile != NULL ) {
+ RefList *rread;
+ rread = read_reflections(rfile);
+ if ( rread == NULL ) {
+ ERROR("Failed to read reference reflections\n");
+ return 1;
+ }
+ reference = asymmetric_indices(rread, sym);
+ reflist_free(rread);
+ ERROR("WARNING: Using an external reference.\n");
+ ERROR("WARNING: If you publish a structure based on the result,"
+ " expect to have to retract your paper!\n");
+ }
+
gsl_set_error_handler_off();
rng = gsl_rng_alloc(gsl_rng_mt19937);
@@ -1008,7 +1203,6 @@ int main(int argc, char *argv[])
n_images = 0;
n_crystals = 0;
n_crystals_seen = 0;
- images = NULL;
crystals = NULL;
if ( sparams_fn != NULL ) {
char line[1024];
@@ -1032,35 +1226,27 @@ int main(int argc, char *argv[])
RefList *as;
int i;
- struct image *images_new;
- struct image *cur;
+ struct image cur;
- images_new = realloc(images, (n_images+1)*sizeof(struct image));
- if ( images_new == NULL ) {
- ERROR("Failed to allocate memory for image list.\n");
- return 1;
- }
- images = images_new;
- cur = &images[n_images];
-
- cur->div = NAN;
- cur->bw = NAN;
- cur->det = NULL;
- if ( read_chunk_2(st, cur, STREAM_READ_REFLECTIONS
- | STREAM_READ_UNITCELL) != 0 ) {
+ cur.div = NAN;
+ cur.bw = NAN;
+ cur.det = NULL;
+ if ( read_chunk_2(st, &cur, STREAM_READ_REFLECTIONS
+ | STREAM_READ_UNITCELL) != 0 ) {
break;
}
- if ( isnan(cur->div) || isnan(cur->bw) ) {
+ if ( isnan(cur.div) || isnan(cur.bw) ) {
ERROR("Chunk doesn't contain beam parameters.\n");
return 1;
}
- for ( i=0; i<cur->n_crystals; i++ ) {
+ for ( i=0; i<cur.n_crystals; i++ ) {
Crystal *cr;
Crystal **crystals_new;
RefList *cr_refl;
+ struct image *image;
n_crystals_seen++;
if ( n_crystals_seen <= start_after ) continue;
@@ -1073,11 +1259,19 @@ int main(int argc, char *argv[])
return 1;
}
crystals = crystals_new;
- crystals[n_crystals] = cur->crystals[i];
+ crystals[n_crystals] = cur.crystals[i];
cr = crystals[n_crystals];
- /* Image pointer will change due to later reallocs */
- crystal_set_image(cr, NULL);
+ image = malloc(sizeof(struct image));
+ if ( image == NULL ) {
+ ERROR("Failed to allocatea memory for image.\n");
+ return 1;
+ }
+
+ crystal_set_image(cr, image);
+ *image = cur;
+ image->n_crystals = 1;
+ image->crystals = &crystals[n_crystals];
/* This is the raw list of reflections */
cr_refl = crystal_get_reflections(cr);
@@ -1087,7 +1281,7 @@ int main(int argc, char *argv[])
if ( polarisation ) {
polarisation_correction(cr_refl,
crystal_get_cell(cr),
- cur);
+ image);
}
if ( !no_free ) select_free_reflections(cr_refl, rng);
@@ -1120,68 +1314,77 @@ int main(int argc, char *argv[])
display_progress(n_images, n_crystals);
fprintf(stderr, "\n");
if ( sparams_fh != NULL ) fclose(sparams_fh);
-
+ audit_info = stream_audit_info(st);
close_stream(st);
- /* Fill in image pointers */
- for ( i=0; i<n_images; i++ ) {
- int j;
- for ( j=0; j<images[i].n_crystals; j++ ) {
-
- Crystal *cryst;
-
- cryst = images[i].crystals[j];
- crystal_set_image(cryst, &images[i]);
-
- /* Now it's safe to do the following */
- update_predictions(cryst);
- calculate_partialities(cryst, pmodel);
-
+ STATUS("Initial partiality calculation...\n");
+ for ( i=0; i<n_crystals; i++ ) {
+ Crystal *cr = crystals[i];
+ if ( force_bandwidth > 0.0 ) {
+ crystal_get_image(cr)->bw = force_bandwidth;
}
+ if ( force_radius > 0.0 ) {
+ crystal_set_profile_radius(cr, force_radius);
+ }
+ update_predictions(cr);
+ calculate_partialities(cr, pmodel);
}
if (csplit != NULL) check_csplit(crystals, n_crystals, csplit);
/* Make a first pass at cutting out crap */
- STATUS("Checking patterns.\n");
+ //STATUS("Early rejection...\n");
//early_rejection(crystals, n_crystals);
- /* Initial rejection, figures of merit etc */
- full = merge_intensities(crystals, n_crystals, nthreads, pmodel,
- min_measurements, push_res, 1);
+ /* Create reference data set if we don't already have one */
+ if ( reference == NULL ) {
+ if ( !no_scale ) {
+ STATUS("Initial scaling...\n");
+ scale_all(crystals, n_crystals, nthreads, scaleflags);
+ }
+ full = merge_intensities(crystals, n_crystals, nthreads,
+ min_measurements, push_res, 1);
+ } else {
+ full = reference;
+ }
+
+ /* Check rejection and write figures of merit */
check_rejection(crystals, n_crystals, full, max_B);
show_all_residuals(crystals, n_crystals, full);
- write_pgraph(full, crystals, n_crystals, 0);
+ if ( !no_pr && !no_logs ) {
+ write_pgraph(full, crystals, n_crystals, 0, "");
+ write_logs_parallel(crystals, n_crystals, full, 0, nthreads,
+ scaleflags);
+ }
/* Iterate */
for ( i=0; i<n_iter; i++ ) {
STATUS("Scaling and refinement cycle %i of %i\n", i+1, n_iter);
- if ( !no_scale ) {
- scale_all(crystals, n_crystals, nthreads, pmodel);
- reflist_free(full);
- full = merge_intensities(crystals, n_crystals, nthreads,
- pmodel, min_measurements,
- push_res, 1);
+ if ( !no_pr ) {
+ refine_all(crystals, n_crystals, full, nthreads, pmodel,
+ 0, i+1, no_logs, sym, amb, scaleflags);
}
- if ( !no_pr ) {
- refine_all(crystals, n_crystals, full, nthreads,
- pmodel);
+ /* Create new reference if needed */
+ if ( reference == NULL ) {
reflist_free(full);
+ if ( !no_scale ) {
+ scale_all(crystals, n_crystals, nthreads,
+ scaleflags);
+ }
full = merge_intensities(crystals, n_crystals, nthreads,
- pmodel, min_measurements,
+ min_measurements,
push_res, 1);
- }
+ } /* else full still equals reference */
check_rejection(crystals, n_crystals, full, max_B);
- reflist_free(full);
- full = merge_intensities(crystals, n_crystals, nthreads,
- pmodel, min_measurements,
- push_res, 1);
show_all_residuals(crystals, n_crystals, full);
- write_pgraph(full, crystals, n_crystals, i+1);
+
+ if ( !no_logs ) {
+ write_pgraph(full, crystals, n_crystals, i+1, "");
+ }
if ( output_everycycle ) {
@@ -1208,21 +1411,37 @@ int main(int argc, char *argv[])
}
}
- /* Dump parameters */
- snprintf(tmp, 1024, "iter%.2d_partialator.params", i+1);
- dump_parameters(tmp, crystals, n_crystals);
+ }
+ }
+ /* Final merge */
+ STATUS("Final merge...\n");
+ if ( reference == NULL ) {
+ reflist_free(full);
+ if ( !no_scale ) {
+ scale_all(crystals, n_crystals, nthreads, scaleflags);
}
+ full = merge_intensities(crystals, n_crystals, nthreads,
+ min_measurements,
+ push_res, 1);
+ } else {
+ full = merge_intensities(crystals, n_crystals, nthreads,
+ min_measurements, push_res, 1);
}
- full = merge_intensities(crystals, n_crystals, nthreads, pmodel,
- min_measurements, push_res, 1);
+ /* Write final figures of merit (no rejection any more) */
show_all_residuals(crystals, n_crystals, full);
- write_pgraph(full, crystals, n_crystals, n_iter+1);
+ if ( !no_pr && !no_logs ) {
+ write_pgraph(full, crystals, n_crystals, -1, "");
+ write_logs_parallel(crystals, n_crystals, full, -1, nthreads,
+ scaleflags);
+ }
/* Output results */
STATUS("Writing overall results to %s\n", outfile);
reflist_add_command_and_version(full, argc, argv);
+ reflist_add_notes(full, "Audit information from stream:");
+ reflist_add_notes(full, audit_info);
write_reflist_2(outfile, full, sym);
/* Output split results */
@@ -1238,9 +1457,6 @@ int main(int argc, char *argv[])
}
}
- /* Dump parameters */
- dump_parameters("partialator.params", crystals, n_crystals);
-
/* Clean up */
gsl_rng_free(rng);
for ( i=0; i<n_crystals; i++ ) {
@@ -1248,13 +1464,9 @@ int main(int argc, char *argv[])
crystal_free(crystals[i]);
}
reflist_free(full);
- free(sym);
+ free_symoplist(sym);
free(outfile);
free(crystals);
- for ( i=0; i<n_images; i++ ) {
- free(images[i].filename);
- }
- free(images);
free(infile);
return 0;
diff --git a/src/post-refinement.c b/src/post-refinement.c
index bf30d299..bde9c02c 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -3,11 +3,11 @@
*
* Post refinement
*
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2015 Thomas White <taw@physics.org>
+ * 2010-2018 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -33,11 +33,8 @@
#include <stdlib.h>
#include <assert.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_vector.h>
-#include <gsl/gsl_linalg.h>
-#include <gsl/gsl_eigen.h>
-#include <gsl/gsl_blas.h>
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_fit.h>
#include "image.h"
#include "post-refinement.h"
@@ -46,15 +43,14 @@
#include "geometry.h"
#include "cell.h"
#include "cell-utils.h"
+#include "reflist-utils.h"
+#include "scaling.h"
+#include "merge.h"
-/* Maximum number of iterations of NLSq to do for each image per macrocycle. */
-#define MAX_CYCLES (10)
-
struct prdata
{
int refined;
- int n_filtered;
};
const char *str_prflag(enum prflag flag)
@@ -79,625 +75,820 @@ const char *str_prflag(enum prflag flag)
case PRFLAG_BIGB :
return "B too big";
+ case PRFLAG_SCALEBAD :
+ return "bad scaling";
+
default :
return "Unknown flag";
}
}
-/* Returns dp(gauss)/dr at "r" */
-static double gaussian_fraction_gradient(double r, double R)
+static UnitCell *rotate_cell_xy(const UnitCell *cell, double ang1, double ang2)
{
- const double ng = 2.6;
- const double sig = R/ng;
-
- /* If the Ewald sphere isn't within the profile, the gradient is zero */
- if ( r < -R ) return 0.0;
- if ( r > +R ) return 0.0;
-
- return exp(-pow(r/sig, 2.0)/2.0) / (sig*sqrt(2.0*M_PI));
+ UnitCell *o;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ double xnew, ynew, znew;
+
+ o = cell_new_from_cell(cell);
+
+ cell_get_reciprocal(o, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+
+ /* "a" around x */
+ xnew = asx;
+ ynew = asy*cos(ang1) + asz*sin(ang1);
+ znew = -asy*sin(ang1) + asz*cos(ang1);
+ asx = xnew; asy = ynew; asz = znew;
+
+ /* "b" around x */
+ xnew = bsx;
+ ynew = bsy*cos(ang1) + bsz*sin(ang1);
+ znew = -bsy*sin(ang1) + bsz*cos(ang1);
+ bsx = xnew; bsy = ynew; bsz = znew;
+
+ /* "c" around x */
+ xnew = csx;
+ ynew = csy*cos(ang1) + csz*sin(ang1);
+ znew = -csy*sin(ang1) + csz*cos(ang1);
+ csx = xnew; csy = ynew; csz = znew;
+
+ /* "a" around y */
+ xnew = asx*cos(ang2) + asz*sin(ang2);
+ ynew = asy;
+ znew = -asx*sin(ang2) + asz*cos(ang2);
+ asx = xnew; asy = ynew; asz = znew;
+
+ /* "b" around y */
+ xnew = bsx*cos(ang2) + bsz*sin(ang2);
+ ynew = bsy;
+ znew = -bsx*sin(ang2) + bsz*cos(ang2);
+ bsx = xnew; bsy = ynew; bsz = znew;
+
+ /* "c" around y */
+ xnew = csx*cos(ang2) + csz*sin(ang2);
+ ynew = csy;
+ znew = -csx*sin(ang2) + csz*cos(ang2);
+ csx = xnew; csy = ynew; csz = znew;
+
+ cell_set_reciprocal(o, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
+
+ return o;
}
-/* Returns dp(sph)/dr at "r" */
-static double sphere_fraction_gradient(double r, double pr)
+static const char *get_label(enum gparam p)
{
- double q, dpdq, dqdr;
-
- /* If the Ewald sphere isn't within the profile, the gradient is zero */
- if ( r < -pr ) return 0.0;
- if ( r > +pr ) return 0.0;
-
- q = (r + pr)/(2.0*pr);
- dpdq = 6.0*(q - q*q);
- dqdr = 1.0 / (2.0*pr);
- return dpdq * dqdr;
+ 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";
+ }
}
-/* Returns dp/dr at "r" */
-static double partiality_gradient(double r, double pr,
- PartialityModel pmodel,
- double rlow, double rhigh)
+/* 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)
{
- double A, D;
-
- D = rlow - rhigh;
-
- switch ( pmodel ) {
+ switch ( p ) {
- default:
- case PMODEL_UNITY:
- return 0.0;
+ 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;
- case PMODEL_SCSPHERE:
- A = sphere_fraction_gradient(r, pr)/D;
- return 4.0*pr*A/3.0;
-
- case PMODEL_SCGAUSSIAN:
- A = gaussian_fraction_gradient(r, pr)/D;
- return 4.0*pr*A/3.0;
+ default : return 0.0;
}
}
-static double sphere_fraction_rgradient(double r, double R)
-{
- /* If the Ewald sphere isn't within the profile, the gradient is zero */
- if ( r < -R ) return 0.0;
- if ( r > +R ) return 0.0;
-
- return -(3.0*r/(4.0*R*R)) * (1.0 - r*r/(R*R));
-}
+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;
-static double gaussian_fraction_rgradient(double r, double R)
-{
- const double ng = 2.6;
- const double sig = R/ng;
+ /* So that it stays around until the end of minimisation */
+ gsl_multimin_function f;
+};
- /* If the Ewald sphere isn't within the profile, the gradient is zero */
- if ( r < -R ) return 0.0;
- if ( r > +R ) return 0.0;
- return -(ng*r/(sqrt(2.0*M_PI)*R*R))*exp(-r*r/(2.0*sig*sig));
+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);
}
-static double volume_fraction_rgradient(double r, double pr,
- PartialityModel pmodel)
+static void apply_parameters(const gsl_vector *v, const gsl_vector *initial,
+ enum gparam *rv, Crystal *cr)
{
- switch ( pmodel )
- {
- case PMODEL_UNITY :
- return 1.0;
-
- case PMODEL_SCSPHERE :
- return sphere_fraction_rgradient(r, pr);
+ int i;
+ double ang1, ang2, R, lambda;
+ UnitCell *cell;
- case PMODEL_SCGAUSSIAN :
- return gaussian_fraction_rgradient(r, pr);
+ /* 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;
- default :
- ERROR("No pmodel in volume_fraction_rgradient!\n");
- return 1.0;
- }
-}
+ for ( i=0; i<v->size; i++ ) {
+ double val;
-static double volume_fraction(double rlow, double rhigh, double pr,
- PartialityModel pmodel)
-{
- switch ( pmodel )
- {
- case PMODEL_UNITY :
- return 1.0;
+ val = get_actual_val(v, initial, rv, i);
- case PMODEL_SCSPHERE :
- return sphere_fraction(rlow, rhigh, pr);
+ switch ( rv[i] ) {
- case PMODEL_SCGAUSSIAN :
- return gaussian_fraction(rlow, rhigh, pr);
+ case GPARAM_ANG1 :
+ ang1 = val;
+ break;
- default :
- ERROR("No pmodel in volume_fraction!\n");
- return 1.0;
- }
-}
+ case GPARAM_ANG2 :
+ ang2 = val;
+ break;
+ case GPARAM_R :
+ R = val;
+ break;
-/* Return the gradient of "fx" wrt parameter 'k' given the current
- * status of the crystal. */
-double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
-{
- double glow, ghigh;
- double rlow, rhigh, p;
- struct image *image = crystal_get_image(cr);
- double R = crystal_get_profile_radius(cr);
- double gr;
+ case GPARAM_WAVELENGTH :
+ lambda = val;
+ break;
- get_partial(refl, &rlow, &rhigh, &p);
+ default :
+ ERROR("Don't understand parameter %i\n", rv[i]);
+ break;
- if ( k == GPARAM_R ) {
+ }
+ }
- double Rglow, Rghigh;
- double D, psph;
+ cell = rotate_cell_xy(crystal_get_cell_const(cr), ang1, ang2);
+ crystal_set_cell(cr, cell);
- D = rlow - rhigh;
- psph = volume_fraction(rlow, rhigh, R, pmodel);
+ crystal_set_profile_radius(cr, R);
+ crystal_get_image(cr)->lambda = lambda;
+}
- Rglow = volume_fraction_rgradient(rlow, R, pmodel);
- Rghigh = volume_fraction_rgradient(rhigh, R, pmodel);
- gr = 4.0*psph/(3.0*D) + (4.0*R/(3.0*D))*(Rglow - Rghigh);
- return gr;
+static double residual_f(const gsl_vector *v, void *pp)
+{
+ struct rf_priv *pv = pp;
+ RefList *list;
+ struct image im;
+ Crystal *cr;
+ UnitCell *cell;
+ double res;
+ int i;
+ for ( i=0; i<v->size; i++ ) {
+ if ( gsl_vector_get(v, i) > 100.0 ) return GSL_NAN;
}
- /* Calculate the gradient of partiality wrt excitation error. */
- glow = partiality_gradient(rlow, R, pmodel, rlow, rhigh);
- ghigh = partiality_gradient(rhigh, R, pmodel, rlow, rhigh);
+ cr = crystal_copy(pv->cr);
+ cell = cell_new_from_cell(crystal_get_cell(cr));
+ if ( cell == NULL ) return GSL_NAN;
+ crystal_set_cell(cr, cell);
+ im = *crystal_get_image(cr);
+ crystal_set_image(cr, &im);
+ apply_parameters(v, pv->initial, pv->rv, cr);
+
+ if ( fabs(crystal_get_profile_radius(cr)) > 5e9 ) {
+ crystal_free(cr);
+ if ( pv->verbose ) STATUS("radius > 5e9\n");
+ return GSL_NAN;
+ }
- if ( k == GPARAM_DIV ) {
+ /* Can happen with grid scans and certain --force-radius values */
+ if ( fabs(crystal_get_profile_radius(cr)) < 0.0000001e9 ) {
+ crystal_free(cr);
+ if ( pv->verbose ) STATUS("radius very small\n");
+ return GSL_NAN;
+ }
- double D, psph, ds;
- signed int hs, ks, ls;
+ if ( im.lambda <= 0.0 ) {
+ crystal_free(cr);
+ if ( pv->verbose ) STATUS("lambda < 0\n");
+ return GSL_NAN;
+ }
- D = rlow - rhigh;
- psph = volume_fraction(rlow, rhigh, R, pmodel);
- get_symmetric_indices(refl, &hs, &ks, &ls);
- ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls);
+ list = copy_reflist(crystal_get_reflections(cr));
+ crystal_set_reflections(cr, list);
- gr = (ds/2.0)*(glow+ghigh) - 4.0*R*psph*ds/(3.0*D*D);
- return gr;
+ update_predictions(cr);
+ calculate_partialities(cr, PMODEL_XSPHERE);
+ if ( scale_one_crystal(cr, pv->full, pv->scaleflags) ) {
+ crystal_free(cr);
+ if ( pv->verbose ) STATUS("Bad scaling\n");
+ return GSL_NAN;
}
+ res = residual(cr, pv->full, 0, NULL, NULL);
+
+ cell_free(crystal_get_cell(cr));
+ reflist_free(crystal_get_reflections(cr));
+ crystal_free(cr);
- gr = r_gradient(crystal_get_cell(cr), k, refl, image) * (glow-ghigh);
- return gr;
+ return res;
}
-static void apply_cell_shift(UnitCell *cell, int k, double shift)
+static double get_initial_param(Crystal *cr, enum gparam p)
{
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
- double as, bs, cs;
+ switch ( p ) {
- cell_get_reciprocal(cell, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz);
- as = modulus(asx, asy, asz);
- bs = modulus(bsx, bsy, bsz);
- cs = modulus(csx, csy, csz);
+ 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;
- /* Forbid any step which looks too large */
- switch ( k )
- {
- case GPARAM_ASX :
- case GPARAM_ASY :
- case GPARAM_ASZ :
- if ( fabs(shift) > 0.1*as ) return;
- break;
-
- case GPARAM_BSX :
- case GPARAM_BSY :
- case GPARAM_BSZ :
- if ( fabs(shift) > 0.1*bs ) return;
- break;
-
- case GPARAM_CSX :
- case GPARAM_CSY :
- case GPARAM_CSZ :
- if ( fabs(shift) > 0.1*cs ) return;
- break;
- }
-
- switch ( k )
- {
- case GPARAM_ASX : asx += shift; break;
- case GPARAM_ASY : asy += shift; break;
- case GPARAM_ASZ : asz += shift; break;
- case GPARAM_BSX : bsx += shift; break;
- case GPARAM_BSY : bsy += shift; break;
- case GPARAM_BSZ : bsz += shift; break;
- case GPARAM_CSX : csx += shift; break;
- case GPARAM_CSY : csy += shift; break;
- case GPARAM_CSZ : csz += shift; break;
- }
-
- cell_set_reciprocal(cell, asx, asy, asz,
- bsx, bsy, bsz,
- csx, csy, csz);
+ default: return 0.0;
+
+ }
}
-/* Apply the given shift to the 'k'th parameter of 'image'. */
-static void apply_shift(Crystal *cr, int k, double shift)
+static int check_angle_shifts(gsl_vector *cur, const gsl_vector *initial,
+ enum gparam *rv, struct rf_priv *residual_f_priv)
{
- double t;
- struct image *image = crystal_get_image(cr);
-
- switch ( k ) {
-
- case GPARAM_DIV :
- if ( shift > 0.1*image->div ) return;
- image->div += shift;
- if ( image->div < 0.0 ) image->div = 0.0;
- break;
-
- case GPARAM_R :
- t = crystal_get_profile_radius(cr);
- if ( shift > 0.1*t ) return;
- t += shift;
- crystal_set_profile_radius(cr, t);
- break;
-
- case GPARAM_ASX :
- case GPARAM_ASY :
- case GPARAM_ASZ :
- case GPARAM_BSX :
- case GPARAM_BSY :
- case GPARAM_BSZ :
- case GPARAM_CSX :
- case GPARAM_CSY :
- case GPARAM_CSZ :
- apply_cell_shift(crystal_get_cell(cr), k, shift);
- break;
+ int i = 0;
+ double ang = 0.0;
- default :
- ERROR("No shift defined for parameter %i\n", k);
- abort();
+ 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;
}
-/* Perform one cycle of post refinement on 'image' against 'full' */
-static double pr_iterate(Crystal *cr, const RefList *full,
- PartialityModel pmodel,
- int *n_filtered, int verbose)
+static RefList *reindex_reflections(RefList *input, SymOpList *amb,
+ SymOpList *sym, int idx)
{
- gsl_matrix *M;
- gsl_vector *v;
- gsl_vector *shifts;
- int param;
+ RefList *n;
Reflection *refl;
RefListIterator *iter;
- RefList *reflections;
- double max_shift;
- int nref = 0;
- int num_params = 0;
- enum gparam rv[32];
- double G, B;
-
- if ( n_filtered != NULL ) *n_filtered = 0;
-
- rv[num_params++] = GPARAM_ASX;
- rv[num_params++] = GPARAM_ASY;
- rv[num_params++] = GPARAM_ASZ;
- rv[num_params++] = GPARAM_BSX;
- rv[num_params++] = GPARAM_BSY;
- rv[num_params++] = GPARAM_BSZ;
- rv[num_params++] = GPARAM_CSX;
- rv[num_params++] = GPARAM_CSY;
- rv[num_params++] = GPARAM_CSZ;
-// rv[num_params++] = GPARAM_R;
-// rv[num_params++] = GPARAM_DIV;
-
- M = gsl_matrix_calloc(num_params, num_params);
- v = gsl_vector_calloc(num_params);
-
- reflections = crystal_get_reflections(cr);
- G = crystal_get_osf(cr);
- B = crystal_get_Bfac(cr);
-
- /* Post-refinement terms */
- for ( refl = first_refl(reflections, &iter);
+
+ n = reflist_new();
+
+ for ( refl = first_refl(input, &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
- signed int ha, ka, la;
- double I_full, delta_I, esd;
- double w;
- double I_partial;
- int k;
- double p, L, s;
- double fx;
- Reflection *match;
- double gradients[num_params];
+ signed int h, k, l;
+ Reflection *rn;
- if ( get_flag(refl) ) continue;
+ get_indices(refl, &h, &k, &l);
+ get_equiv(amb, NULL, idx, h, k, l, &h, &k, &l);
+ get_asymm(sym, h, k, l, &h, &k, &l);
+ rn = add_refl(n, h, k, l);
+ copy_data(rn, refl);
+
+ get_symmetric_indices(rn, &h, &k, &l);
+ get_equiv(amb, NULL, idx, h, k, l, &h, &k, &l);
+ set_symmetric_indices(rn, h, k, l);
+ }
- get_indices(refl, &ha, &ka, &la);
- match = find_refl(full, ha, ka, la);
- if ( match == NULL ) continue;
- I_full = get_intensity(match);
+ return n;
+}
- if ( get_redundancy(match) < 2 ) continue;
- p = get_partiality(refl);
- L = get_lorentz(refl);
- I_partial = get_intensity(refl);
- esd = get_esd_intensity(refl);
- s = resolution(crystal_get_cell(cr), ha, ka, la);
+static void reindex_cell(UnitCell *cell, SymOpList *amb, int idx)
+{
+ signed int h, k, l;
+ struct rvec na, nb, nc;
+ struct rvec as, bs, cs;
+
+ cell_get_reciprocal(cell, &as.u, &as.v, &as.w,
+ &bs.u, &bs.v, &bs.w,
+ &cs.u, &cs.v, &cs.w);
+
+ get_equiv(amb, NULL, idx, 1, 0, 0, &h, &k, &l);
+ na.u = as.u*h + bs.u*k + cs.u*l;
+ na.v = as.v*h + bs.v*k + cs.v*l;
+ na.w = as.w*h + bs.w*k + cs.w*l;
+
+ get_equiv(amb, NULL, idx, 0, 1, 0, &h, &k, &l);
+ nb.u = as.u*h + bs.u*k + cs.u*l;
+ nb.v = as.v*h + bs.v*k + cs.v*l;
+ nb.w = as.w*h + bs.w*k + cs.w*l;
+
+ get_equiv(amb, NULL, idx, 0, 0, 1, &h, &k, &l);
+ nc.u = as.u*h + bs.u*k + cs.u*l;
+ nc.v = as.v*h + bs.v*k + cs.v*l;
+ nc.w = as.w*h + bs.w*k + cs.w*l;
+
+ cell_set_reciprocal(cell, na.u, na.v, na.w,
+ nb.u, nb.v, nb.w,
+ nc.u, nc.v, nc.w);
+}
- if ( I_partial < 3.0*esd ) continue;
- /* Calculate the weight for this reflection */
- w = (s/1e9)*(s/1e9) / (esd*esd);
+static void try_reindex(Crystal *crin, const RefList *full,
+ SymOpList *sym, SymOpList *amb, int scaleflags)
+{
+ RefList *list;
+ Crystal *cr;
+ UnitCell *cell;
+ double residual_original, residual_flipped;
+ int idx, n;
- /* Calculate all gradients for this reflection */
- for ( k=0; k<num_params; k++ ) {
- gradients[k] = gradient(cr, rv[k], refl, pmodel);
- gradients[k] *= exp(G)*exp(-B*s*s)*I_full/L;
- }
+ if ( sym == NULL || amb == NULL ) return;
- for ( k=0; k<num_params; k++ ) {
+ if ( scale_one_crystal(crin, full, scaleflags) ) return;
+ residual_original = residual(crin, full, 0, NULL, NULL);
- int g;
- double v_c, v_curr;
+ cr = crystal_copy(crin);
- for ( g=0; g<num_params; g++ ) {
+ n = num_equivs(amb, NULL);
- double M_c, M_curr;
+ for ( idx=0; idx<n; idx++ ) {
- /* Matrix is symmetric */
- if ( g > k ) continue;
+ cell = cell_new_from_cell(crystal_get_cell(crin));
+ if ( cell == NULL ) return;
+ reindex_cell(cell, amb, idx);
+ crystal_set_cell(cr, cell);
- M_c = w * gradients[g] * gradients[k];
+ list = reindex_reflections(crystal_get_reflections(crin),
+ amb, sym, idx);
+ crystal_set_reflections(cr, list);
- M_curr = gsl_matrix_get(M, k, g);
- gsl_matrix_set(M, k, g, M_curr + M_c);
- gsl_matrix_set(M, g, k, M_curr + M_c);
+ update_predictions(cr);
+ calculate_partialities(cr, PMODEL_XSPHERE);
+
+ if ( scale_one_crystal(cr, full, scaleflags) ) return;
+ residual_flipped = residual(cr, full, 0, NULL, NULL);
+
+ if ( residual_flipped < residual_original ) {
+ crystal_set_cell(crin, cell);
+ crystal_set_reflections(crin, list);
+ residual_original = residual_flipped;
+ } else {
+ cell_free(crystal_get_cell(cr));
+ reflist_free(crystal_get_reflections(cr));
+ }
+ }
- }
+ crystal_free(cr);
+}
- fx = exp(G)*p*exp(-B*s*s)*I_full/L;
- delta_I = I_partial - fx;
- v_c = w * delta_I * gradients[k];
- v_curr = gsl_vector_get(v, k);
- gsl_vector_set(v, k, v_curr + v_c);
- }
+void write_test_logs(Crystal *crystal, const RefList *full,
+ signed int cycle, int serial)
+{
+ FILE *fh;
+ struct image *image = crystal_get_image(crystal);
+ char tmp[256];
+ char ins[5];
- nref++;
- }
- if ( verbose ) {
- //STATUS("The original equation:\n");
- //show_matrix_eqn(M, v);
- STATUS("%i reflections went into the equations.\n", nref);
- }
-
- if ( nref < num_params ) {
- crystal_set_user_flag(cr, PRFLAG_FEWREFL);
- gsl_matrix_free(M);
- gsl_vector_free(v);
- return 0.0;
- }
-
- max_shift = 0.0;
- shifts = solve_svd(v, M, n_filtered, 0);
- if ( shifts != NULL ) {
-
- for ( param=0; param<num_params; param++ ) {
- double shift = gsl_vector_get(shifts, param);
- if ( verbose ) STATUS("Shift %i: %e\n", param, shift);
- if ( isnan(shift) ) {
- //ERROR("NaN shift parameter %i (image ser %i), "
- // "%i reflections.\n", rv[param],
- // crystal_get_image(cr)->serial,
- // nref);
- } else {
- apply_shift(cr, rv[param], shift);
- }
- if ( fabs(shift) > max_shift ) max_shift = fabs(shift);
- }
+ snprintf(tmp, 256, "pr-logs/parameters-crystal%i.dat", serial);
+ if ( cycle == 0 ) {
+ fh = fopen(tmp, "w");
} else {
- crystal_set_user_flag(cr, PRFLAG_SOLVEFAIL);
+ fh = fopen(tmp, "a");
}
- gsl_matrix_free(M);
- gsl_vector_free(v);
- gsl_vector_free(shifts);
+ if ( fh == NULL ) {
+ ERROR("Failed to open '%s'\n", tmp);
+ return;
+ }
+
+ if ( cycle == 0 ) {
+ fprintf(fh, "Image: %s %s\n",
+ image->filename, get_event_string(image->event));
+ }
+
+ if ( cycle >= 0 ) {
+ snprintf(ins, 4, "%i", cycle);
+ } else {
+ ins[0] = 'F';
+ ins[1] = '\0';
+ }
- return max_shift;
+ fprintf(fh, "%s rlp_size = %e\n", ins, crystal_get_profile_radius(crystal)/1e10);
+ fprintf(fh, "%s mosaicity = %e\n", ins, crystal_get_mosaicity(crystal));
+ fprintf(fh, "%s wavelength = %f A\n", ins, image->lambda*1e10);
+ fprintf(fh, "%s bandwidth = %f\n", ins, image->bw);
+ fprintf(fh, "%s my scale factor = %e\n", ins, crystal_get_osf(crystal));
+
+ double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
+ cell_get_reciprocal(crystal_get_cell(crystal), &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+ fprintf(fh, "%s %f, %f, %f, %f, %f, %f, %f, %f, %f\n", ins,
+ asx/1e10, bsx/1e10, csx/1e10,
+ asy/1e10, bsy/1e10, csy/1e10,
+ asz/1e10, bsz/1e10, csz/1e10);
+ fclose(fh);
}
-double residual(Crystal *cr, const RefList *full, int free,
- int *pn_used, const char *filename)
+void write_specgraph(Crystal *crystal, const RefList *full,
+ signed int cycle, int serial)
{
- double dev = 0.0;
- double G, B;
+ FILE *fh;
+ char tmp[256];
Reflection *refl;
RefListIterator *iter;
- FILE *fh = NULL;
- int n_used = 0;
+ double G = crystal_get_osf(crystal);
+ double B = crystal_get_Bfac(crystal);
+ UnitCell *cell;
+ struct image *image = crystal_get_image(crystal);
+ char ins[5];
- if ( filename != NULL ) {
- fh = fopen(filename, "a");
- if ( fh == NULL ) {
- ERROR("Failed to open '%s'\n", filename);
- }
+ snprintf(tmp, 256, "pr-logs/specgraph-crystal%i.dat", serial);
+
+ if ( cycle == 0 ) {
+ fh = fopen(tmp, "w");
+ } else {
+ fh = fopen(tmp, "a");
+ }
+
+ if ( fh == NULL ) {
+ ERROR("Failed to open '%s'\n", tmp);
+ return;
+ }
+
+ if ( cycle == 0 ) {
+ fprintf(fh, "Image: %s %s\n",
+ image->filename, get_event_string(image->event));
+ fprintf(fh, "khalf/m 1/d(m) pcalc pobs iteration h k l\n");
}
- G = crystal_get_osf(cr);
- B = crystal_get_Bfac(cr);
+ cell = crystal_get_cell(crystal);
- for ( refl = first_refl(crystal_get_reflections(cr), &iter);
+ if ( cycle >= 0 ) {
+ snprintf(ins, 4, "%i", cycle);
+ } else {
+ ins[0] = 'F';
+ ins[1] = '\0';
+ }
+
+ for ( refl = first_refl(crystal_get_reflections(crystal), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
- double p, L, s, w;
+ double corr, Ipart, Ifull, pobs, pcalc;
+ double res;
signed int h, k, l;
Reflection *match;
- double esd, I_full, I_partial;
- double fx, dc;
-
- if ( free && !get_flag(refl) ) continue;
get_indices(refl, &h, &k, &l);
+ res = resolution(cell, h, k, l);
+
match = find_refl(full, h, k, l);
if ( match == NULL ) continue;
- I_full = get_intensity(match);
- if ( get_redundancy(match) < 2 ) continue;
+ /* Don't calculate pobs if reference reflection is weak */
+ if ( fabs(get_intensity(match)) / get_esd_intensity(match) < 3.0 ) continue;
- p = get_partiality(refl);
- L = get_lorentz(refl);
- I_partial = get_intensity(refl);
- esd = get_esd_intensity(refl);
- s = resolution(crystal_get_cell(cr), h, k, l);
+ corr = G * exp(B*res*res) * get_lorentz(refl);
+ Ipart = get_intensity(refl) * corr;
+ Ifull = get_intensity(match);
+ pobs = Ipart / Ifull;
+ pcalc = get_partiality(refl);
- if ( I_partial < 3.0*esd ) continue;
+ fprintf(fh, "%e %e %f %f %s %4i %4i %4i\n",
+ get_khalf(refl), 2.0*res, pcalc, pobs, ins, h, k, l);
- fx = exp(G)*p*exp(-B*s*s)*I_full/L;
- dc = I_partial - fx;
- w = (s/1e9)*(s/1e9)/(esd*esd);
- dev += w*dc*dc;
- n_used++;
-
- if ( fh != NULL ) {
- fprintf(fh, "%4i %4i %4i %e %e\n",
- h, k, l, s, dev);
- }
}
- if ( fh != NULL ) fclose(fh);
-
- if ( pn_used != NULL ) *pn_used = n_used;
- return dev;
+ fclose(fh);
}
-static void write_residual_graph(Crystal *cr, const RefList *full)
+static gsl_multimin_fminimizer *setup_minimiser(Crystal *cr, const RefList *full,
+ int verbose, int serial,
+ int scaleflags,
+ struct rf_priv *priv)
{
+ gsl_multimin_fminimizer *min;
+ int n_params = 0;
int i;
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
- double a;
- double step;
- double orig_asx;
- FILE *fh;
- UnitCell *cell;
- cell = crystal_get_cell(cr);
+ /* 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; i<n_params; i++ ) {
+ gsl_vector_set(priv->initial, i, get_initial_param(cr, priv->rv[i]));
+ gsl_vector_set(priv->vals, i, 0.0);
+ gsl_vector_set(priv->step, i, 1.0);
+ }
- fh = fopen("residual-plot.dat", "w");
+ min = gsl_multimin_fminimizer_alloc(gsl_multimin_fminimizer_nmsimplex2,
+ n_params);
+ gsl_multimin_fminimizer_set(min, &priv->f, priv->vals, priv->step);
- cell_get_reciprocal(cell, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz);
- a = modulus(asx, asy, asz);
- orig_asx = asx;
+ return min;
+}
- step = a/100.0;
- for ( i=-50; i<=50; i++ ) {
+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)
+{
+ FILE *fh;
+ char fn[64];
+ char ins[5];
+ gsl_multimin_fminimizer *min;
+ struct rf_priv priv;
+ int idx1, idx2;
+ int i;
- double res;
- int n;
- asx = orig_asx + (i*step);
- cell_set_reciprocal(cell, asx, asy, asz,
- bsx, bsy, bsz,
- csx, csy, csz);
- update_predictions(cr);
- calculate_partialities(cr, PMODEL_SCSPHERE);
- res = residual(cr, full, 0, &n, NULL);
- fprintf(fh, "%i %e %e %i\n", i, asx, res, n);
+ min = setup_minimiser(cr, full, 0, serial, scaleflags, &priv);
+
+ idx1 = 99;
+ idx2 = 99;
+ for ( i=0; i<priv.f.n; i++ ) {
+ if ( priv.rv[i] == par1 ) idx1 = i;
+ if ( priv.rv[i] == par2 ) idx2 = i;
}
+ assert(idx1 != 99);
+ assert(idx2 != 99);
- cell_set_reciprocal(cell, orig_asx, asy, asz,
- bsx, bsy, bsz,
- csx, csy, csz);
- update_predictions(cr);
- calculate_partialities(cr, PMODEL_SCSPHERE);
+ if ( cycle >= 0 ) {
+ snprintf(ins, 4, "%i", cycle);
+ } else {
+ ins[0] = 'F';
+ ins[1] = '\0';
+ }
+
+ snprintf(fn, 63, "pr-logs/grid-crystal%i-cycle%s-%s.dat",
+ serial, ins, name);
+ 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 ) {
+ int first=1;
+ for ( v1=-5.0; v1<=5.0; v1+=0.25 ) {
+ double res;
+ gsl_vector_set(min->x, idx1, v1);
+ gsl_vector_set(min->x, idx2, v2);
+ res = residual_f(min->x, &priv);
+ if ( !first ) fprintf(fh, " ");
+ first = 0;
+ fprintf(fh, "%e", res);
+ }
+ fprintf(fh, "\n");
+ }
+ }
fclose(fh);
+
+ gsl_multimin_fminimizer_free(min);
+ gsl_vector_free(priv.initial);
+ gsl_vector_free(priv.vals);
+ gsl_vector_free(priv.step);
+}
+
+
+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");
}
static void do_pr_refine(Crystal *cr, const RefList *full,
- PartialityModel pmodel, int verbose)
+ PartialityModel pmodel, int verbose, int serial,
+ int cycle, int write_logs,
+ SymOpList *sym, SymOpList *amb, int scaleflags)
{
- int i, done;
- double old_dev;
- UnitCell *cell = crystal_get_cell(cr);
+ gsl_multimin_fminimizer *min;
+ struct rf_priv priv;
+ int n_iter = 0;
+ int status;
+ double residual_start, residual_free_start;
+ FILE *fh = NULL;
+
+ try_reindex(cr, full, sym, amb, scaleflags);
+
+ if ( scale_one_crystal(cr, full, scaleflags | SCALE_VERBOSE_ERRORS) ) {
+ 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);
+ }
- old_dev = residual(cr, full, 0, 0, NULL);
+ min = setup_minimiser(cr, full, verbose, serial, scaleflags, &priv);
if ( verbose ) {
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
- cell_get_reciprocal(cell, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz);
- STATUS("Initial asx = %e\n", asx);
- STATUS("PR initial dev = %10.5e, free dev = %10.5e\n",
- old_dev, residual(cr, full, 1, NULL, NULL));
- }
-
- i = 0;
- done = 0;
- do {
+ 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);
+ }
- double dev;
+ if ( write_logs ) {
- pr_iterate(cr, full, pmodel, NULL, verbose);
+ char fn[64];
- update_predictions(cr);
- calculate_partialities(cr, pmodel);
+ 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");
+ 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);
+ }
- dev = residual(cr, full, 0, 0, NULL);
- if ( fabs(dev - old_dev) < dev*0.0001 ) done = 1;
+ }
+
+ do {
+ double res;
+
+ n_iter++;
+
+ status = gsl_multimin_fminimizer_iterate(min);
+ if ( status ) break;
+
+ res = residual_f(min->x, &priv);
+ if ( isnan(res) ) {
+ status = GSL_ENOPROG;
+ break;
+ }
if ( verbose ) {
- STATUS("PR iter %2i: dev = %10.5e, free dev = %10.5e\n",
- i+1, dev, residual(cr, full, 1, NULL, NULL));
+ 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);
}
- i++;
- old_dev = dev;
+ 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);
+ }
+
+ status = gsl_multimin_test_size(min->size, 0.005);
- } while ( i < 30 && !done );
+ } while ( status == GSL_CONTINUE && n_iter < 1000 );
if ( verbose ) {
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
- cell_get_reciprocal(cell, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz);
- STATUS("Final asx = %e\n", asx);
+ STATUS("Done with refinement after %i iter\n", n_iter);
+ STATUS("status = %i (%s)\n", status, gsl_strerror(status));
}
-}
+ if ( status == GSL_SUCCESS ) {
-static struct prdata pr_refine(Crystal *cr, const RefList *full,
- PartialityModel pmodel)
-{
- int verbose = 0;
- struct prdata prdata;
+ if ( check_angle_shifts(min->x, priv.initial, priv.rv, &priv) ) return;
- prdata.refined = 0;
- prdata.n_filtered = 0;
+ if ( verbose ) {
- if ( verbose ) {
- write_residual_graph(cr, full);
+ 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);
+ fclose(fh);
+ }
+
+ /* Apply the final shifts */
+ apply_parameters(min->x, priv.initial, priv.rv, cr);
+ 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\n", serial);
}
- do_pr_refine(cr, full, pmodel, verbose);
+ if ( write_logs ) {
+ write_gridscan(cr, full, cycle, serial, scaleflags);
+ write_specgraph(cr, full, cycle, serial);
+ write_test_logs(cr, full, cycle, serial);
+ }
- if ( crystal_get_user_flag(cr) == 0 ) {
- prdata.refined = 1;
+ if ( crystal_get_profile_radius(cr) > 5e9 ) {
+ ERROR("Very large radius: crystal %i\n", serial);
}
- return prdata;
+ gsl_multimin_fminimizer_free(min);
+ gsl_vector_free(priv.initial);
+ gsl_vector_free(priv.vals);
+ gsl_vector_free(priv.step);
}
@@ -706,7 +897,14 @@ struct refine_args
RefList *full;
Crystal *crystal;
PartialityModel pmodel;
+ int serial;
struct prdata prdata;
+ int verbose;
+ int cycle;
+ int no_logs;
+ SymOpList *sym;
+ SymOpList *amb;
+ int scaleflags;
};
@@ -724,8 +922,18 @@ static void refine_image(void *task, int id)
{
struct refine_args *pargs = task;
Crystal *cr = pargs->crystal;
+ 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,
+ pargs->serial, pargs->cycle, write_logs,
+ pargs->sym, pargs->amb, pargs->scaleflags);
- pargs->prdata = pr_refine(cr, pargs->full, pargs->pmodel);
+ if ( crystal_get_user_flag(cr) == 0 ) {
+ pargs->prdata.refined = 1;
+ }
}
@@ -738,6 +946,7 @@ static void *get_image(void *vqargs)
memcpy(task, &qargs->task_defaults, sizeof(struct refine_args));
task->crystal = qargs->crystals[qargs->n_started];
+ task->serial = qargs->n_started;
qargs->n_started++;
@@ -757,7 +966,9 @@ static void done_image(void *vqargs, void *task)
void refine_all(Crystal **crystals, int n_crystals,
- RefList *full, int nthreads, PartialityModel pmodel)
+ RefList *full, int nthreads, PartialityModel pmodel,
+ int verbose, int cycle, int no_logs,
+ SymOpList *sym, SymOpList *amb, int scaleflags)
{
struct refine_args task_defaults;
struct queue_args qargs;
@@ -766,7 +977,12 @@ void refine_all(Crystal **crystals, int n_crystals,
task_defaults.crystal = NULL;
task_defaults.pmodel = pmodel;
task_defaults.prdata.refined = 0;
- task_defaults.prdata.n_filtered = 0;
+ task_defaults.verbose = verbose;
+ task_defaults.cycle = cycle;
+ task_defaults.no_logs = no_logs;
+ task_defaults.sym = sym;
+ task_defaults.amb = amb;
+ task_defaults.scaleflags = scaleflags;
qargs.task_defaults = task_defaults;
qargs.n_started = 0;
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 7f395013..71a6d7f3 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -3,11 +3,11 @@
*
* Post refinement
*
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2015 Thomas White <taw@physics.org>
+ * 2010-2018 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -41,6 +41,7 @@
#include "utils.h"
#include "crystal.h"
#include "geometry.h"
+#include "symmetry.h"
enum prflag
@@ -51,19 +52,28 @@ enum prflag
PRFLAG_EARLY = 18,
PRFLAG_CC = 19,
PRFLAG_BIGB = 20,
+ PRFLAG_SCALEBAD = 21,
};
extern const char *str_prflag(enum prflag flag);
extern void refine_all(Crystal **crystals, int n_crystals,
- RefList *full, int nthreads, PartialityModel pmodel);
+ RefList *full, int nthreads, PartialityModel pmodel,
+ int verbose, int cycle, int no_logs,
+ SymOpList *sym, SymOpList *amb, int scaleflags);
+
+extern void write_gridscan(Crystal *cr, const RefList *full,
+ int cycle, int serial, int scaleflags);
+
+extern void write_specgraph(Crystal *crystal, const RefList *full,
+ signed int cycle, int serial);
/* Exported so it can be poked by tests/pr_p_gradient_check */
extern double gradient(Crystal *cr, int k, Reflection *refl,
PartialityModel pmodel);
-extern double residual(Crystal *cr, const RefList *full, int free,
- int *pn_used, const char *filename);
+extern void write_test_logs(Crystal *crystal, const RefList *full,
+ signed int cycle, int serial);
#endif /* POST_REFINEMENT_H */
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 5930c189..ab5b9af1 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -519,6 +519,7 @@ int main(int argc, char *argv[])
double push_res = +INFINITY;
double min_cc = -INFINITY;
int twopass = 0;
+ char *audit_info;
/* Long options */
const struct option longopts[] = {
@@ -822,11 +823,14 @@ int main(int argc, char *argv[])
hist_nbins);
}
+ audit_info = stream_audit_info(st);
+ close_stream(st);
+
reflist_add_command_and_version(model, argc, argv);
+ reflist_add_notes(model, "Audit information from stream:");
+ reflist_add_notes(model, audit_info);
write_reflist_2(output, model, sym);
- close_stream(st);
-
free_symoplist(sym);
reflist_free(model);
free(output);
diff --git a/src/process_image.c b/src/process_image.c
index 87089289..4b02e694 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -301,7 +301,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
/* Integrate! */
time_accounts_set(taccs, TACC_INTEGRATION);
sb_shared->pings[cookie]++;
- integrate_all_5(&image, iargs->int_meth, PMODEL_SCSPHERE,
+ integrate_all_5(&image, iargs->int_meth, PMODEL_XSPHERE,
iargs->push_res,
iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
iargs->int_diag, iargs->int_diag_h,
diff --git a/src/rejection.c b/src/rejection.c
index 10f13f9c..9d41f9b4 100644
--- a/src/rejection.c
+++ b/src/rejection.c
@@ -137,7 +137,7 @@ static void check_cc(Crystal *cr, RefList *full)
pcalc = get_partiality(refl);
/* Observed partiality */
- corr = exp(-G) * exp(B*res*res) * get_lorentz(refl);
+ corr = G * exp(B*res*res) * get_lorentz(refl);
Ipart = get_intensity(refl) * corr;
pobs = Ipart / get_intensity(match);
diff --git a/src/scaling.c b/src/scaling.c
index 77a64384..cabc5952 100644
--- a/src/scaling.c
+++ b/src/scaling.c
@@ -3,11 +3,11 @@
*
* Scaling
*
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2015 Thomas White <taw@physics.org>
+ * 2010-2017 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -38,290 +38,22 @@
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_blas.h>
+#include <gsl/gsl_fit.h>
+#include <gsl/gsl_statistics_double.h>
#include "merge.h"
#include "post-refinement.h"
#include "symmetry.h"
#include "cell.h"
#include "cell-utils.h"
-
-
-/* Maximum number of iterations of NLSq to do for each image per macrocycle. */
-#define MAX_CYCLES (10)
-
-
-/* Apply the given shift to the 'k'th parameter of 'image'. */
-static void apply_shift(Crystal *cr, int k, double shift)
-{
- double t;
-
- switch ( k ) {
-
- case GPARAM_BFAC :
- t = crystal_get_Bfac(cr);
- t += shift;
- crystal_set_Bfac(cr, t);
- break;
-
- case GPARAM_OSF :
- t = crystal_get_osf(cr);
- t += shift;
- crystal_set_osf(cr, t);
- break;
-
- default :
- ERROR("No shift defined for parameter %i\n", k);
- abort();
-
- }
-}
-
-
-/* Perform one cycle of scaling of 'cr' against 'full' */
-static double scale_iterate(Crystal *cr, const RefList *full,
- PartialityModel pmodel, int *nr)
-{
- gsl_matrix *M;
- gsl_vector *v;
- gsl_vector *shifts;
- int param;
- Reflection *refl;
- RefListIterator *iter;
- RefList *reflections;
- double max_shift;
- int nref = 0;
- int num_params = 0;
- enum gparam rv[32];
- double G, B;
-
- *nr = 0;
-
- rv[num_params++] = GPARAM_OSF;
- rv[num_params++] = GPARAM_BFAC;
-
- M = gsl_matrix_calloc(num_params, num_params);
- v = gsl_vector_calloc(num_params);
-
- reflections = crystal_get_reflections(cr);
- G = crystal_get_osf(cr);
- B = crystal_get_Bfac(cr);
-
- /* Scaling terms */
- for ( refl = first_refl(reflections, &iter);
- refl != NULL;
- refl = next_refl(refl, iter) )
- {
- signed int ha, ka, la;
- double I_full, delta_I, esd;
- double w;
- double I_partial;
- int k;
- double p, L, s;
- double fx;
- Reflection *match;
- double gradients[num_params];
-
- /* If reflection is free-flagged, don't use it here */
- if ( get_flag(refl) ) continue;
-
- /* Find the full version */
- get_indices(refl, &ha, &ka, &la);
- match = find_refl(full, ha, ka, la);
- if ( match == NULL ) continue;
-
- /* Merged intensitty */
- I_full = get_intensity(match);
-
- /* Actual measurement of this reflection from this pattern */
- I_partial = get_intensity(refl);
- esd = get_esd_intensity(refl);
- p = get_partiality(refl);
-
- /* Scale only using strong reflections */
- if ( I_partial <= 3.0*esd ) continue; /* Also because of log */
- if ( get_redundancy(match) < 2 ) continue;
- if ( I_full <= 0 ) continue; /* Because log */
- if ( p <= 0.0 ) continue; /* Because of log */
-
- L = get_lorentz(refl);
- s = resolution(crystal_get_cell(cr), ha, ka, la);
-
- /* Calculate the weight for this reflection */
- w = 1.0;
-
- /* Calculate all gradients for this reflection */
- for ( k=0; k<num_params; k++ ) {
-
- if ( rv[k] == GPARAM_OSF ) {
- gradients[k] = 1.0;
- } else if ( rv[k] == GPARAM_BFAC ) {
- gradients[k] = -s*s;
- } else {
- ERROR("Unrecognised scaling gradient.\n");
- abort();
- }
- }
-
- for ( k=0; k<num_params; k++ ) {
-
- int g;
- double v_c, v_curr;
-
- for ( g=0; g<num_params; g++ ) {
-
- double M_c, M_curr;
-
- /* Matrix is symmetric */
- if ( g > k ) continue;
-
- M_c = w * gradients[g] * gradients[k];
-
- M_curr = gsl_matrix_get(M, k, g);
- gsl_matrix_set(M, k, g, M_curr + M_c);
- gsl_matrix_set(M, g, k, M_curr + M_c);
-
- }
-
- fx = G + log(p) - log(L) - B*s*s + log(I_full);
- delta_I = log(I_partial) - fx;
- v_c = w * delta_I * gradients[k];
- v_curr = gsl_vector_get(v, k);
- gsl_vector_set(v, k, v_curr + v_c);
-
- }
-
- nref++;
- }
-
- *nr = nref;
-
- if ( nref < num_params ) {
- crystal_set_user_flag(cr, PRFLAG_FEWREFL);
- gsl_matrix_free(M);
- gsl_vector_free(v);
- return 0.0;
- }
-
- max_shift = 0.0;
- shifts = solve_svd(v, M, NULL, 0);
- if ( shifts != NULL ) {
-
- for ( param=0; param<num_params; param++ ) {
- double shift = gsl_vector_get(shifts, param);
- apply_shift(cr, rv[param], shift);
- if ( fabs(shift) > max_shift ) max_shift = fabs(shift);
- }
-
- } else {
- crystal_set_user_flag(cr, PRFLAG_SOLVEFAIL);
- }
-
- gsl_matrix_free(M);
- gsl_vector_free(v);
- gsl_vector_free(shifts);
-
- return max_shift;
-}
-
-
-double log_residual(Crystal *cr, const RefList *full, int free,
- int *pn_used, const char *filename)
-{
- double dev = 0.0;
- double G, B;
- Reflection *refl;
- RefListIterator *iter;
- int n_used = 0;
- FILE *fh = NULL;
-
- G = crystal_get_osf(cr);
- B = crystal_get_Bfac(cr);
- if ( filename != NULL ) {
- fh = fopen(filename, "a");
- if ( fh == NULL ) {
- ERROR("Failed to open '%s'\n", filename);
- }
- }
-
- for ( refl = first_refl(crystal_get_reflections(cr), &iter);
- refl != NULL;
- refl = next_refl(refl, iter) )
- {
- double p, L, s, w;
- signed int h, k, l;
- Reflection *match;
- double esd, I_full, I_partial;
- double fx, dc;
-
- if ( free && !get_flag(refl) ) continue;
-
- get_indices(refl, &h, &k, &l);
- match = find_refl(full, h, k, l);
- if ( match == NULL ) continue;
-
- p = get_partiality(refl);
- L = get_lorentz(refl);
- I_partial = get_intensity(refl);
- I_full = get_intensity(match);
- esd = get_esd_intensity(refl);
- s = resolution(crystal_get_cell(cr), h, k, l);
-
- if ( I_partial <= 3.0*esd ) continue; /* Also because of log */
- if ( get_redundancy(match) < 2 ) continue;
- if ( I_full <= 0 ) continue; /* Because log */
- if ( p <= 0.0 ) continue; /* Because of log */
-
- fx = G + log(p) - log(L) - B*s*s + log(I_full);
- dc = log(I_partial) - fx;
- w = 1.0;
- dev += w*dc*dc;
-
- if ( fh != NULL ) {
- fprintf(fh, "%4i %4i %4i %e %e\n",
- h, k, l, s, dev);
- }
-
- }
-
- if ( fh != NULL ) fclose(fh);
-
- if ( pn_used != NULL ) *pn_used = n_used;
- return dev;
-}
-
-
-static void do_scale_refine(Crystal *cr, const RefList *full,
- PartialityModel pmodel, int *nr)
-{
- int i, done;
- double old_dev;
-
- old_dev = log_residual(cr, full, 0, NULL, NULL);
-
- i = 0;
- done = 0;
- do {
-
- double dev;
-
- scale_iterate(cr, full, pmodel, nr);
-
- dev = log_residual(cr, full, 0, 0, NULL);
- if ( fabs(dev - old_dev) < dev*0.01 ) done = 1;
-
- i++;
- old_dev = dev;
-
- } while ( i < MAX_CYCLES && !done );
-}
+#include "scaling.h"
struct scale_args
{
RefList *full;
Crystal *crystal;
- PartialityModel pmodel;
- int n_reflections;
+ int flags;
};
@@ -331,7 +63,6 @@ struct queue_args
int n_done;
Crystal **crystals;
int n_crystals;
- long long int n_reflections;
struct scale_args task_defaults;
};
@@ -339,8 +70,7 @@ struct queue_args
static void scale_crystal(void *task, int id)
{
struct scale_args *pargs = task;
- do_scale_refine(pargs->crystal, pargs->full, pargs->pmodel,
- &pargs->n_reflections);
+ scale_one_crystal(pargs->crystal, pargs->full, pargs->flags);
}
@@ -363,11 +93,7 @@ static void *get_crystal(void *vqargs)
static void done_crystal(void *vqargs, void *task)
{
struct queue_args *qa = vqargs;
- struct scale_args *wargs = task;
-
qa->n_done++;
- qa->n_reflections += wargs->n_reflections;
-
progress_bar(qa->n_done, qa->n_crystals, "Scaling");
free(task);
}
@@ -381,14 +107,12 @@ static double total_log_r(Crystal **crystals, int n_crystals, RefList *full,
int n = 0;
for ( i=0; i<n_crystals; i++ ) {
-
double r;
if ( crystal_get_user_flag(crystals[i]) ) continue;
r = log_residual(crystals[i], full, 0, NULL, NULL);
if ( isnan(r) ) continue;
total += r;
n++;
-
}
if ( ninc != NULL ) *ninc = n;
@@ -397,8 +121,7 @@ static double total_log_r(Crystal **crystals, int n_crystals, RefList *full,
/* Perform iterative scaling, all the way to convergence */
-void scale_all(Crystal **crystals, int n_crystals, int nthreads,
- PartialityModel pmodel)
+void scale_all(Crystal **crystals, int n_crystals, int nthreads, int scaleflags)
{
struct scale_args task_defaults;
struct queue_args qargs;
@@ -406,7 +129,7 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads,
int niter = 0;
task_defaults.crystal = NULL;
- task_defaults.pmodel = pmodel;
+ task_defaults.flags = scaleflags;
qargs.task_defaults = task_defaults;
qargs.n_crystals = n_crystals;
@@ -422,18 +145,15 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads,
double bef_res;
full = merge_intensities(crystals, n_crystals, nthreads,
- pmodel, 2, INFINITY, 0);
+ 2, INFINITY, 0);
old_res = new_res;
bef_res = total_log_r(crystals, n_crystals, full, NULL);
qargs.task_defaults.full = full;
qargs.n_started = 0;
qargs.n_done = 0;
- qargs.n_reflections = 0;
run_threads(nthreads, scale_crystal, get_crystal, done_crystal,
&qargs, n_crystals, 0, 0, 0);
- STATUS("%lli reflections went into the scaling.\n",
- qargs.n_reflections);
new_res = total_log_r(crystals, n_crystals, full, &ninc);
STATUS("Log residual went from %e to %e, %i crystals\n",
@@ -456,3 +176,164 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads,
ERROR("Too many iterations - giving up!\n");
}
}
+
+
+/* Calculates G and B, by which cr's reflections should be multiplied to fit reference */
+int scale_one_crystal(Crystal *cr, const RefList *listR, int flags)
+{
+ const Reflection *reflS;
+ RefListIterator *iter;
+ int max_n = 256;
+ int n = 0;
+ double *x;
+ double *y;
+ double *w;
+ int r;
+ double cov00, cov01, cov11, chisq;
+ int n_esdS = 0;
+ int n_esdR = 0;
+ int n_ihS = 0;
+ int n_ihR = 0;
+ int n_nanS = 0;
+ int n_nanR = 0;
+ int n_infS = 0;
+ int n_infR = 0;
+ int n_part = 0;
+ int n_nom = 0;
+ int n_red = 0;
+ RefList *listS = crystal_get_reflections(cr);
+ UnitCell *cell = crystal_get_cell(cr);
+ double G, B;
+
+ assert(cell != NULL);
+ assert(listR != NULL);
+ assert(listS != NULL);
+
+ x = malloc(max_n*sizeof(double));
+ w = malloc(max_n*sizeof(double));
+ y = malloc(max_n*sizeof(double));
+ if ( (x==NULL) || (y==NULL) || (w==NULL) ) {
+ ERROR("Failed to allocate memory for scaling.\n");
+ return 1;
+ }
+
+ int nb = 0;
+ for ( reflS = first_refl_const(listS, &iter);
+ reflS != NULL;
+ reflS = next_refl_const(reflS, iter) )
+ {
+ signed int h, k, l;
+ const Reflection *reflR;
+ double IhR, IhS, esdS, pS, LS;
+ double s;
+
+ nb++;
+
+ get_indices(reflS, &h, &k, &l);
+ reflR = find_refl(listR, h, k, l);
+ if ( reflR == NULL ) {
+ n_nom++;
+ continue;
+ }
+
+ s = resolution(cell, h, k, l);
+
+ IhR = get_intensity(reflR);
+ IhS = get_intensity(reflS);
+ esdS = get_esd_intensity(reflS);
+ pS = get_partiality(reflS);
+ LS = get_lorentz(reflS);
+
+ /* Problem cases in approximate descending order of severity */
+ if ( isnan(IhR) ) { n_nanR++; continue; }
+ if ( isinf(IhR) ) { n_infR++; continue; }
+ if ( isnan(IhS) ) { n_nanS++; continue; }
+ if ( isinf(IhS) ) { n_infS++; continue; }
+ if ( pS < 0.3 ) { n_part++; continue; }
+ if ( IhS <= 0.0 ) { n_ihS++; continue; }
+ if ( IhS <= 3.0*esdS ) { n_esdS++; continue; }
+ if ( IhR <= 0.0 ) { n_ihR++; continue; }
+ if ( get_redundancy(reflR) < 2 ) { n_red++; continue; }
+
+ if ( n == max_n ) {
+ max_n *= 2;
+ x = realloc(x, max_n*sizeof(double));
+ y = realloc(y, max_n*sizeof(double));
+ w = realloc(w, max_n*sizeof(double));
+ if ( (x==NULL) || (y==NULL) || (w==NULL) ) {
+ ERROR("Failed to allocate memory for scaling.\n");
+ return 1;
+ }
+ }
+
+ x[n] = s*s;
+ y[n] = log(LS) + log(IhS) -log(pS) - log(IhR);
+ w[n] = pS;
+ n++;
+
+ }
+
+ if ( n < 2 ) {
+ if ( flags & SCALE_VERBOSE_ERRORS ) {
+ ERROR("Not enough reflections for scaling (had %i, but %i remain)\n", nb, n);
+ if ( n_esdR ) ERROR("%i reference reflection esd\n", n_esdR);
+ if ( n_esdS ) ERROR("%i subject reflection esd\n", n_esdS);
+ if ( n_ihR ) ERROR("%i reference reflection intensity\n", n_ihR);
+ if ( n_red ) ERROR("%i reference reflection redundancy\n", n_red);
+ if ( n_ihS ) ERROR("%i subject reflection intensity\n", n_ihS);
+ if ( n_nanR ) ERROR("%i reference reflection nan\n", n_nanR);
+ if ( n_nanS ) ERROR("%i subject reflection nan\n", n_nanS);
+ if ( n_infR ) ERROR("%i reference reflection inf\n", n_infR);
+ if ( n_infS ) ERROR("%i subject reflection inf\n", n_infS);
+ if ( n_part ) ERROR("%i subject reflection partiality\n", n_part);
+ if ( n_nom ) ERROR("%i no match in reference list\n", n_nom);
+ }
+ free(x);
+ free(y);
+ free(w);
+ return 1;
+ }
+
+ if ( flags & SCALE_NO_B ) {
+ G = gsl_stats_wmean(w, 1, y, 1, n);
+ B = 0.0;
+ r = 0;
+ } else {
+ r = gsl_fit_wlinear(x, 1, w, 1, y, 1, n, &G, &B, &cov00, &cov01, &cov11, &chisq);
+ }
+
+ if ( r ) {
+ ERROR("Scaling failed.\n");
+ free(x);
+ free(y);
+ free(w);
+ return 1;
+ }
+
+ if ( isnan(G) ) {
+
+ if ( flags & SCALE_VERBOSE_ERRORS ) {
+ ERROR("Scaling gave NaN (%i pairs)\n", n);
+ if ( n < 10 ) {
+ int i;
+ for ( i=0; i<n; i++ ) {
+ STATUS("%3i %e %e %e\n", i, x[i], y[i], w[i]);
+ }
+ }
+ }
+
+ free(x);
+ free(y);
+ free(w);
+ return 1;
+ }
+
+ crystal_set_osf(cr, exp(G));
+ crystal_set_Bfac(cr, -B);
+
+ free(x);
+ free(y);
+ free(w);
+
+ return 0;
+}
diff --git a/src/scaling.h b/src/scaling.h
index 534043d9..8c2e7ef2 100644
--- a/src/scaling.h
+++ b/src/scaling.h
@@ -3,11 +3,11 @@
*
* Scaling
*
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2015 Thomas White <taw@physics.org>
+ * 2010-2018 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -38,10 +38,16 @@
#include "crystal.h"
#include "geometry.h"
-extern double log_residual(Crystal *cr, const RefList *full, int free,
- int *pn_used, const char *filename);
+enum ScaleFlags
+{
+ SCALE_NONE = 0,
+ SCALE_NO_B = 1<<0, /* Don't use Debye-Waller part */
+ SCALE_VERBOSE_ERRORS = 1<<1,
+};
+
+extern int scale_one_crystal(Crystal *cr, const RefList *reference, int flags);
extern void scale_all(Crystal **crystals, int n_crystals, int nthreads,
- PartialityModel pmodel);
+ int flags);
#endif /* SCALING_H */
generated by cgit v1.2.3 (git 2.25.1) at 2024-08-13 03:27:42 +0000