From bc628cd85f512d962828e9b8c5df5a6e6d22fac5 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 7 Mar 2013 10:50:13 +0100 Subject: Add the possibility to have different partiality models --- libcrystfel/src/geometry.c | 2 +- libcrystfel/src/geometry.h | 18 +++++++++++++++++- src/partialator.c | 25 ++++++++++++++++++++++--- src/post-refinement.c | 13 +++++++------ src/post-refinement.h | 6 ++++-- tests/pr_gradient_check.c | 25 +++++++++++++------------ 6 files changed, 64 insertions(+), 25 deletions(-) diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index 621c9484..450ebfaf 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -346,7 +346,7 @@ RefList *select_intersections(struct image *image, Crystal *cryst) /* Calculate partialities and apply them to the image's reflections */ -void update_partialities(Crystal *cryst) +void update_partialities(Crystal *cryst, PartialityModel pmodel) { Reflection *refl; RefListIterator *iter; diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h index aecdc28a..8b726ed7 100644 --- a/libcrystfel/src/geometry.h +++ b/libcrystfel/src/geometry.h @@ -43,10 +43,26 @@ extern "C" { #endif +/** + * PartialityModel: + * @PMODEL_SPHERE : Intersection of sphere with excited volume of reciprocal + * space. + * @PMODEL_UNITY : Set all all partialities and Lorentz factors to 1. + * + * A %PartialityModelModel describes a geometrical model which can be used to + * calculate spot partialities and Lorentz correction factors. + **/ +typedef enum { + + PMODEL_SPHERE, + PMODEL_UNITY, + +} PartialityModel; + extern RefList *find_intersections(struct image *image, Crystal *cryst); extern RefList *select_intersections(struct image *image, Crystal *cryst); -extern void update_partialities(Crystal *cryst); +extern void update_partialities(Crystal *cryst, PartialityModel pmodel); #ifdef __cplusplus } diff --git a/src/partialator.c b/src/partialator.c index a4af3c18..079e078e 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -88,6 +88,7 @@ struct refine_args { RefList *full; Crystal *crystal; + PartialityModel pmodel; }; @@ -106,7 +107,7 @@ static void refine_image(void *task, int id) struct refine_args *pargs = task; Crystal *cr = pargs->crystal; - pr_refine(cr, pargs->full); + pr_refine(cr, pargs->full, pargs->pmodel); } @@ -314,6 +315,8 @@ int main(int argc, char *argv[]) int noscale = 0; Stream *st; Crystal **crystals; + char *pmodel_str = NULL; + PartialityModel pmodel = PMODEL_SPHERE; /* Long options */ const struct option longopts[] = { @@ -326,6 +329,7 @@ int main(int argc, char *argv[]) {"iterations", 1, NULL, 'n'}, {"no-scale", 0, &noscale, 1}, {"reference", 1, NULL, 'r'}, + {"partiality", 1, NULL, 'm'}, {0, 0, NULL, 0} }; @@ -336,7 +340,7 @@ int main(int argc, char *argv[]) } /* Short options */ - while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:r:j:", + while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:r:j:m:", longopts, NULL)) != -1) { @@ -370,6 +374,10 @@ int main(int argc, char *argv[]) n_iter = atoi(optarg); break; + case 'm' : + pmodel_str = strdup(optarg); + break; + case 'b' : beam = get_beam_parameters(optarg); if ( beam == NULL ) { @@ -434,6 +442,17 @@ int main(int argc, char *argv[]) return 1; } + if ( pmodel_str != NULL ) { + if ( strcmp(pmodel_str, "sphere") == 0 ) { + pmodel = PMODEL_SPHERE; + } else if ( strcmp(pmodel_str, "unity") == 0 ) { + pmodel = PMODEL_UNITY; + } else { + ERROR("Unknown partiality model '%s'.\n", pmodel_str); + return 1; + } + } + if ( reference_file != NULL ) { RefList *list; @@ -544,7 +563,7 @@ int main(int argc, char *argv[]) crystal_set_image(cryst, &images[i]); /* Now it's safe to do the following */ - update_partialities(cryst); + update_partialities(cryst, pmodel); as = crystal_get_reflections(cryst); nobs += select_scalable_reflections(as, reference); diff --git a/src/post-refinement.c b/src/post-refinement.c index 1439b148..f6f4caa5 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -88,7 +88,7 @@ static double partiality_rgradient(double r, double profile_radius) /* Return the gradient of parameter 'k' given the current status of 'image'. */ -double gradient(Crystal *cr, int k, Reflection *refl) +double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) { double ds, azix, aziy; double ttlow, tthigh, tt; @@ -364,7 +364,8 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M) /* Perform one cycle of post refinement on 'image' against 'full' */ -static double pr_iterate(Crystal *cr, const RefList *full) +static double pr_iterate(Crystal *cr, const RefList *full, + PartialityModel pmodel) { gsl_matrix *M; gsl_vector *v; @@ -420,7 +421,7 @@ static double pr_iterate(Crystal *cr, const RefList *full) /* Calculate all gradients for this reflection */ for ( k=0; k Date: Thu, 7 Mar 2013 11:16:23 +0100 Subject: Add PartialityModel to docs --- doc/reference/libcrystfel/CrystFEL-sections.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/doc/reference/libcrystfel/CrystFEL-sections.txt b/doc/reference/libcrystfel/CrystFEL-sections.txt index 720cabff..ff8001cf 100644 --- a/doc/reference/libcrystfel/CrystFEL-sections.txt +++ b/doc/reference/libcrystfel/CrystFEL-sections.txt @@ -380,6 +380,7 @@ crystal_set_user_flag
geometry +PartialityModel find_intersections select_intersections update_partialities -- cgit v1.2.3 From aedd1501af8707b9f502a2fd765d743f53c3d078 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 7 Mar 2013 13:23:41 +0100 Subject: pr_gradient_check: Fix and improve --- tests/pr_gradient_check.c | 72 +++++++++++++++++++++++++++++------------------ 1 file changed, 44 insertions(+), 28 deletions(-) diff --git a/tests/pr_gradient_check.c b/tests/pr_gradient_check.c index c7eb0129..2631b1bb 100644 --- a/tests/pr_gradient_check.c +++ b/tests/pr_gradient_check.c @@ -155,6 +155,7 @@ static Crystal *new_shifted_crystal(Crystal *cr, int refine, double incr_val) case REF_R : cell = cell_new_from_cell(crystal_get_cell(cr)); + crystal_set_cell(cr_new, cell); crystal_set_profile_radius(cr_new, r + incr_val); break; @@ -218,8 +219,8 @@ static void calc_either_side(Crystal *cr, double incr_val, -static int test_gradients(Crystal *cr, double incr_val, int refine, - const char *str, PartialityModel pmodel) +static double test_gradients(Crystal *cr, double incr_val, int refine, + const char *str, PartialityModel pmodel) { Reflection *refl; RefListIterator *iter; @@ -227,9 +228,11 @@ static int test_gradients(Crystal *cr, double incr_val, int refine, int i; int *valid; int nref; - int n_acc, n_valid; + int n_good, n_invalid, n_small, n_nan, n_bad; RefList *reflections; //FILE *fh; + int ntot = 0; + double total = 0.0; reflections = find_intersections(crystal_get_image(cr), cr); crystal_set_reflections(cr, reflections); @@ -261,7 +264,8 @@ static int test_gradients(Crystal *cr, double incr_val, int refine, //fh = fopen("wrongness.dat", "a"); - n_valid = nref; n_acc = 0; + n_invalid = 0; n_good = 0; + n_nan = 0; n_small = 0; n_bad = 0; i = 0; for ( refl = first_refl(reflections, &iter); refl != NULL; @@ -275,7 +279,7 @@ static int test_gradients(Crystal *cr, double incr_val, int refine, get_indices(refl, &h, &k, &l); if ( !valid[i] ) { - n_valid--; + n_invalid++; } else { double r1, r2, p; @@ -289,8 +293,20 @@ static int test_gradients(Crystal *cr, double incr_val, int refine, get_partial(refl, &r1, &r2, &p, &cl, &ch); - if ( (fabs(cgrad) > 5e-8) && - !within_tolerance(grad, cgrad, 10.0) ) + if ( fabs(cgrad) < 5e-8 ) { + n_small++; + continue; + } + + if ( isnan(cgrad) ) { + n_nan++; + continue; + } + + total += fabs(cgrad - grad); + ntot++; + + if ( !within_tolerance(grad, cgrad, 5.0) ) { STATUS("!- %s %3i %3i %3i" @@ -298,6 +314,7 @@ static int test_gradients(Crystal *cr, double incr_val, int refine, " %10.2e %10.2e\n", str, h, k, l, grad, cgrad, cgrad/grad, r1, r2); + n_bad++; } else { @@ -307,7 +324,7 @@ static int test_gradients(Crystal *cr, double incr_val, int refine, // str, h, k, l, grad, cgrad, cgrad/grad, // r1, r2); - n_acc++; + n_good++; } @@ -323,13 +340,12 @@ static int test_gradients(Crystal *cr, double incr_val, int refine, } - STATUS("%s: %i out of %i valid gradients were accurate.\n", - str, n_acc, n_valid); + STATUS("%3s: %3i within 5%%, %3i outside, %3i nan, %3i invalid, " + "%3i small. ", str, n_good, n_bad, n_nan, n_invalid, n_small); //fclose(fh); - if ( n_acc != n_valid ) return 1; - - return 0; + STATUS("Fractional error = %f %%\n", 100.0*total/ntot); + return total / ntot; } @@ -345,8 +361,8 @@ int main(int argc, char *argv[]) Crystal *cr; struct quaternion orientation; int i; - int val; const PartialityModel pmodel = PMODEL_SPHERE; + double tot = 0.0; image.width = 1024; image.height = 1024; @@ -374,8 +390,6 @@ int main(int argc, char *argv[]) deg2rad(90.0), deg2rad(90.0)); - val = 0; - for ( i=0; i<1; i++ ) { UnitCell *rot; @@ -389,32 +403,34 @@ int main(int argc, char *argv[]) &bz, &cx, &cy, &cz); incr_val = incr_frac * image.div; - val += test_gradients(cr, incr_val, REF_DIV, "div", pmodel); + tot += test_gradients(cr, incr_val, REF_DIV, "div", pmodel); + + incr_val = incr_frac * crystal_get_profile_radius(cr); + tot += test_gradients(cr, incr_val, REF_R, "R", pmodel); incr_val = incr_frac * ax; - val += test_gradients(cr, incr_val, REF_ASX, "ax*", pmodel); + tot += test_gradients(cr, incr_val, REF_ASX, "ax*", pmodel); incr_val = incr_frac * ay; - val += test_gradients(cr, incr_val, REF_ASY, "ay*", pmodel); + tot += test_gradients(cr, incr_val, REF_ASY, "ay*", pmodel); incr_val = incr_frac * az; - val += test_gradients(cr, incr_val, REF_ASZ, "az*", pmodel); + tot += test_gradients(cr, incr_val, REF_ASZ, "az*", pmodel); incr_val = incr_frac * bx; - val += test_gradients(cr, incr_val, REF_BSX, "bx*", pmodel); + tot += test_gradients(cr, incr_val, REF_BSX, "bx*", pmodel); incr_val = incr_frac * by; - val += test_gradients(cr, incr_val, REF_BSY, "by*", pmodel); + tot += test_gradients(cr, incr_val, REF_BSY, "by*", pmodel); incr_val = incr_frac * bz; - val += test_gradients(cr, incr_val, REF_BSZ, "bz*", pmodel); + tot += test_gradients(cr, incr_val, REF_BSZ, "bz*", pmodel); incr_val = incr_frac * cx; - val += test_gradients(cr, incr_val, REF_CSX, "cx*", pmodel); + tot += test_gradients(cr, incr_val, REF_CSX, "cx*", pmodel); incr_val = incr_frac * cy; - val += test_gradients(cr, incr_val, REF_CSY, "cy*", pmodel); + tot += test_gradients(cr, incr_val, REF_CSY, "cy*", pmodel); incr_val = incr_frac * cz; - val += test_gradients(cr, incr_val, REF_CSZ, "cz*", pmodel); + tot += test_gradients(cr, incr_val, REF_CSZ, "cz*", pmodel); } - STATUS("Returning 0 by default: gradients only needed for experimental" - " features of CrystFEL.\n"); + STATUS("Mean fractional error = %f %%\n", 100.0*tot/10.0); return 0; } -- cgit v1.2.3 From f8de3ad3f3410b95180f569d99b9b1f9bd9523ab Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 8 Mar 2013 10:31:29 +0100 Subject: pr_gradient_check: Look for correlations, not individual values --- tests/pr_gradient_check.c | 103 ++++++++++++++++++++++++++++------------------ 1 file changed, 64 insertions(+), 39 deletions(-) diff --git a/tests/pr_gradient_check.c b/tests/pr_gradient_check.c index 2631b1bb..cf11d18c 100644 --- a/tests/pr_gradient_check.c +++ b/tests/pr_gradient_check.c @@ -29,6 +29,7 @@ #include #include +#include #include #include @@ -230,9 +231,14 @@ static double test_gradients(Crystal *cr, double incr_val, int refine, int nref; int n_good, n_invalid, n_small, n_nan, n_bad; RefList *reflections; - //FILE *fh; + FILE *fh; int ntot = 0; double total = 0.0; + char tmp[32]; + double *vec1; + double *vec2; + int n_line; + double cc; reflections = find_intersections(crystal_get_image(cr), cr); crystal_set_reflections(cr, reflections); @@ -240,7 +246,7 @@ static double test_gradients(Crystal *cr, double incr_val, int refine, nref = num_reflections(reflections); if ( nref < 10 ) { ERROR("Too few reflections found. Failing test by default.\n"); - return -1; + return 0.0; } vals[0] = malloc(nref*sizeof(long double)); @@ -248,13 +254,13 @@ static double test_gradients(Crystal *cr, double incr_val, int refine, vals[2] = malloc(nref*sizeof(long double)); if ( (vals[0] == NULL) || (vals[1] == NULL) || (vals[2] == NULL) ) { ERROR("Couldn't allocate memory.\n"); - return -1; + return 0.0; } valid = malloc(nref*sizeof(int)); if ( valid == NULL ) { ERROR("Couldn't allocate memory.\n"); - return -1; + return 0.0; } for ( i=0; iindexed_cell, h, k, l), - //rad2deg(tt), - // cgrad, - // fabs((grad-cgrad)/grad)); - } i++; @@ -342,10 +355,11 @@ static double test_gradients(Crystal *cr, double incr_val, int refine, STATUS("%3s: %3i within 5%%, %3i outside, %3i nan, %3i invalid, " "%3i small. ", str, n_good, n_bad, n_nan, n_invalid, n_small); - //fclose(fh); + fclose(fh); - STATUS("Fractional error = %f %%\n", 100.0*total/ntot); - return total / ntot; + cc = gsl_stats_correlation(vec1, 1, vec2, 1, n_line); + STATUS("CC = %+f\n", cc); + return cc; } @@ -362,7 +376,7 @@ int main(int argc, char *argv[]) struct quaternion orientation; int i; const PartialityModel pmodel = PMODEL_SPHERE; - double tot = 0.0; + int fail = 0; image.width = 1024; image.height = 1024; @@ -393,6 +407,7 @@ int main(int argc, char *argv[]) for ( i=0; i<1; i++ ) { UnitCell *rot; + double val; orientation = random_quaternion(); rot = cell_rotate(cell, orientation); @@ -403,34 +418,44 @@ int main(int argc, char *argv[]) &bz, &cx, &cy, &cz); incr_val = incr_frac * image.div; - tot += test_gradients(cr, incr_val, REF_DIV, "div", pmodel); + val = test_gradients(cr, incr_val, REF_DIV, "div", pmodel); + if ( val > 0.1 ) fail = 1; incr_val = incr_frac * crystal_get_profile_radius(cr); - tot += test_gradients(cr, incr_val, REF_R, "R", pmodel); + val = test_gradients(cr, incr_val, REF_R, "R", pmodel); + if ( val > 0.1 ) fail = 1; incr_val = incr_frac * ax; - tot += test_gradients(cr, incr_val, REF_ASX, "ax*", pmodel); + val = test_gradients(cr, incr_val, REF_ASX, "ax*", pmodel); + if ( val > 0.1 ) fail = 1; incr_val = incr_frac * ay; - tot += test_gradients(cr, incr_val, REF_ASY, "ay*", pmodel); + val = test_gradients(cr, incr_val, REF_ASY, "ay*", pmodel); + if ( val > 0.1 ) fail = 1; incr_val = incr_frac * az; - tot += test_gradients(cr, incr_val, REF_ASZ, "az*", pmodel); + val = test_gradients(cr, incr_val, REF_ASZ, "az*", pmodel); + if ( val > 0.1 ) fail = 1; incr_val = incr_frac * bx; - tot += test_gradients(cr, incr_val, REF_BSX, "bx*", pmodel); + val = test_gradients(cr, incr_val, REF_BSX, "bx*", pmodel); + if ( val > 0.1 ) fail = 1; incr_val = incr_frac * by; - tot += test_gradients(cr, incr_val, REF_BSY, "by*", pmodel); + val = test_gradients(cr, incr_val, REF_BSY, "by*", pmodel); + if ( val > 0.1 ) fail = 1; incr_val = incr_frac * bz; - tot += test_gradients(cr, incr_val, REF_BSZ, "bz*", pmodel); + val = test_gradients(cr, incr_val, REF_BSZ, "bz*", pmodel); + if ( val > 0.1 ) fail = 1; incr_val = incr_frac * cx; - tot += test_gradients(cr, incr_val, REF_CSX, "cx*", pmodel); + val = test_gradients(cr, incr_val, REF_CSX, "cx*", pmodel); + if ( val > 0.1 ) fail = 1; incr_val = incr_frac * cy; - tot += test_gradients(cr, incr_val, REF_CSY, "cy*", pmodel); + val = test_gradients(cr, incr_val, REF_CSY, "cy*", pmodel); + if ( val > 0.1 ) fail = 1; incr_val = incr_frac * cz; - tot += test_gradients(cr, incr_val, REF_CSZ, "cz*", pmodel); + val = test_gradients(cr, incr_val, REF_CSZ, "cz*", pmodel); + if ( val > 0.1 ) fail = 1; } - STATUS("Mean fractional error = %f %%\n", 100.0*tot/10.0); - return 0; + return fail; } -- cgit v1.2.3 From 3a348d3f7e2586440590747c1b921b8ce6dbc0e1 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 8 Mar 2013 10:49:49 +0100 Subject: Work on post refinement Brought across from "tom/pr" Conflicts: src/indexamajig.c src/post-refinement.c tests/pr_gradient_check.c --- libcrystfel/src/utils.h | 20 ++++++++++++++ src/post-refinement.c | 69 +++++++++++++++++++++++++---------------------- tests/pr_gradient_check.c | 3 +++ 3 files changed, 60 insertions(+), 32 deletions(-) diff --git a/libcrystfel/src/utils.h b/libcrystfel/src/utils.h index b75693db..1adb69e6 100644 --- a/libcrystfel/src/utils.h +++ b/libcrystfel/src/utils.h @@ -139,6 +139,11 @@ static inline double modulus(double x, double y, double z) return sqrt(x*x + y*y + z*z); } +static inline double modulus2d(double x, double y) +{ + return sqrt(x*x + y*y); +} + static inline double modulus_squared(double x, double y, double z) { return x*x + y*y + z*z; } @@ -165,6 +170,21 @@ static inline double angle_between(double x1, double y1, double z1, return acos(cosine); } +/* Answer in radians */ +static inline double angle_between_2d(double x1, double y1, + double x2, double y2) +{ + double mod1 = modulus2d(x1, y1); + double mod2 = modulus2d(x2, y2); + double cosine = (x1*x2 + y1*y2) / (mod1*mod2); + + /* Fix domain if outside due to rounding */ + if ( cosine > 1.0 ) cosine = 1.0; + if ( cosine < -1.0 ) cosine = -1.0; + + return acos(cosine); +} + static inline int within_tolerance(double a, double b, double percent) { double tol = fabs(a) * (percent/100.0); diff --git a/src/post-refinement.c b/src/post-refinement.c index f6f4caa5..9e63736f 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -90,8 +90,7 @@ static double partiality_rgradient(double r, double profile_radius) /* Return the gradient of parameter 'k' given the current status of 'image'. */ double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) { - double ds, azix, aziy; - double ttlow, tthigh, tt; + double ds, azi; double nom, den; double g; double asx, asy, asz; @@ -99,9 +98,11 @@ double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) double csx, csy, csz; double xl, yl, zl; signed int hs, ks, ls; - double r1, r2, p; + double rlow, rhigh, p; int clamp_low, clamp_high; - double klow, khigh; + double philow, phihigh, phi; + double khigh, klow; + double tl, cet, cez; double gr; struct image *image = crystal_get_image(cr); double r = crystal_get_profile_radius(cr); @@ -116,33 +117,37 @@ double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) zl = hs*asz + ks*bsz + ls*csz; ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls); - get_partial(refl, &r1, &r2, &p, &clamp_low, &clamp_high); + get_partial(refl, &rlow, &rhigh, &p, &clamp_low, &clamp_high); + /* "low" gives the largest Ewald sphere (wavelength short => k large) + * "high" gives the smallest Ewald sphere (wavelength long => k small) + */ klow = 1.0/(image->lambda - image->lambda*image->bw/2.0); khigh = 1.0/(image->lambda + image->lambda*image->bw/2.0); - ttlow = angle_between(0.0, 0.0, 1.0, xl, yl, zl+klow); - tthigh = angle_between(0.0, 0.0, 1.0, xl, yl, zl+khigh); - if ( (clamp_low == 0) && (clamp_high == 0) ) { - tt = (ttlow+tthigh)/2.0; - } else if ( clamp_high == 0 ) { - tt = tthigh + image->div; - } else if ( clamp_low == 0 ) { - tt = ttlow - image->div; - } else { - tt = 0.0; - /* Gradient should come out as zero in this case */ - } - azix = angle_between(1.0, 0.0, 0.0, xl, yl, 0.0); - aziy = angle_between(0.0, 1.0, 0.0, xl, yl, 0.0); + tl = sqrt(xl*xl + yl*yl); + ds = modulus(xl, yl, zl); + + cet = -sin(image->div/2.0) * klow; + cez = -cos(image->div/2.0) * klow; + philow = M_PI_2 - angle_between_2d(tl-cet, zl-cez, 1.0, 0.0); + + cet = -sin(image->div/2.0) * khigh; + cez = -cos(image->div/2.0) * khigh; + phihigh = M_PI_2 - angle_between_2d(tl-cet, zl-cez, 1.0, 0.0); + + /* Approximation: philow and phihigh are very similar */ + phi = (philow + phihigh) / 2.0; + + azi = atan2(yl, xl); /* Calculate the gradient of partiality wrt excitation error. */ g = 0.0; if ( clamp_low == 0 ) { - g -= partiality_gradient(r1, r); + g -= partiality_gradient(rlow, r); } if ( clamp_high == 0 ) { - g += partiality_gradient(r2, r); + g += partiality_gradient(rhigh, r); } /* For many gradients, just multiply the above number by the gradient @@ -167,40 +172,40 @@ double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) case REF_R : g = 0.0; if ( clamp_low == 0 ) { - g += partiality_rgradient(r1, r); + g += partiality_rgradient(rlow, r); } if ( clamp_high == 0 ) { - g += partiality_rgradient(r2, r); + g += partiality_rgradient(rhigh, r); } return g; /* Cell parameters and orientation */ case REF_ASX : - return hs * sin(tt) * cos(azix) * g; + return hs * sin(phi) * cos(azi) * g; case REF_BSX : - return ks * sin(tt) * cos(azix) * g; + return ks * sin(phi) * cos(azi) * g; case REF_CSX : - return ls * sin(tt) * cos(azix) * g; + return ls * sin(phi) * cos(azi) * g; case REF_ASY : - return hs * sin(tt) * cos(aziy) * g; + return hs * sin(phi) * sin(azi) * g; case REF_BSY : - return ks * sin(tt) * cos(aziy) * g; + return ks * sin(phi) * sin(azi) * g; case REF_CSY : - return ls * sin(tt) * cos(aziy) * g; + return ls * sin(phi) * sin(azi) * g; case REF_ASZ : - return hs * cos(tt) * g; + return hs * cos(phi) * g; case REF_BSZ : - return ks * cos(tt) * g; + return ks * cos(phi) * g; case REF_CSZ : - return ls * cos(tt) * g; + return ls * cos(phi) * g; } diff --git a/tests/pr_gradient_check.c b/tests/pr_gradient_check.c index cf11d18c..d5344801 100644 --- a/tests/pr_gradient_check.c +++ b/tests/pr_gradient_check.c @@ -425,6 +425,9 @@ int main(int argc, char *argv[]) val = test_gradients(cr, incr_val, REF_R, "R", pmodel); if ( val > 0.1 ) fail = 1; + //incr_val = incr_frac * image.profile_radius; + //val += test_gradients(&image, incr_val, REF_R, "rad"); + incr_val = incr_frac * ax; val = test_gradients(cr, incr_val, REF_ASX, "ax*", pmodel); if ( val > 0.1 ) fail = 1; -- cgit v1.2.3 From d268cf95f3fc8de4b349c429e5c24448f0e85bcb Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 8 Mar 2013 10:53:46 +0100 Subject: Group x, y and z gradients together --- tests/pr_gradient_check.c | 40 +++++++++++++++++++++------------------- 1 file changed, 21 insertions(+), 19 deletions(-) diff --git a/tests/pr_gradient_check.c b/tests/pr_gradient_check.c index d5344801..0dd5466b 100644 --- a/tests/pr_gradient_check.c +++ b/tests/pr_gradient_check.c @@ -221,7 +221,8 @@ static void calc_either_side(Crystal *cr, double incr_val, static double test_gradients(Crystal *cr, double incr_val, int refine, - const char *str, PartialityModel pmodel) + const char *str, const char *file, + PartialityModel pmodel) { Reflection *refl; RefListIterator *iter; @@ -268,7 +269,7 @@ static double test_gradients(Crystal *cr, double incr_val, int refine, calc_either_side(cr, incr_val, valid, vals, refine); - snprintf(tmp, 32, "gradient-test-%s.dat", str); + snprintf(tmp, 32, "gradient-test-%s.dat", file); fh = fopen(tmp, "w"); vec1 = malloc(nref*sizeof(double)); @@ -418,44 +419,45 @@ int main(int argc, char *argv[]) &bz, &cx, &cy, &cz); incr_val = incr_frac * image.div; - val = test_gradients(cr, incr_val, REF_DIV, "div", pmodel); + val = test_gradients(cr, incr_val, REF_DIV, "div", "div", + pmodel); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * crystal_get_profile_radius(cr); - val = test_gradients(cr, incr_val, REF_R, "R", pmodel); + val = test_gradients(cr, incr_val, REF_R, "R", "R", pmodel); if ( val > 0.1 ) fail = 1; //incr_val = incr_frac * image.profile_radius; //val += test_gradients(&image, incr_val, REF_R, "rad"); incr_val = incr_frac * ax; - val = test_gradients(cr, incr_val, REF_ASX, "ax*", pmodel); + val = test_gradients(cr, incr_val, REF_ASX, "ax*", "x", pmodel); if ( val > 0.1 ) fail = 1; - incr_val = incr_frac * ay; - val = test_gradients(cr, incr_val, REF_ASY, "ay*", pmodel); + incr_val = incr_frac * bx; + val = test_gradients(cr, incr_val, REF_BSX, "bx*", "x", pmodel); if ( val > 0.1 ) fail = 1; - incr_val = incr_frac * az; - val = test_gradients(cr, incr_val, REF_ASZ, "az*", pmodel); + incr_val = incr_frac * cx; + val = test_gradients(cr, incr_val, REF_CSX, "cx*", "x", pmodel); if ( val > 0.1 ) fail = 1; - incr_val = incr_frac * bx; - val = test_gradients(cr, incr_val, REF_BSX, "bx*", pmodel); + incr_val = incr_frac * ay; + val = test_gradients(cr, incr_val, REF_ASY, "ay*", "y", pmodel); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * by; - val = test_gradients(cr, incr_val, REF_BSY, "by*", pmodel); + val = test_gradients(cr, incr_val, REF_BSY, "by*", "y", pmodel); if ( val > 0.1 ) fail = 1; - incr_val = incr_frac * bz; - val = test_gradients(cr, incr_val, REF_BSZ, "bz*", pmodel); + incr_val = incr_frac * cy; + val = test_gradients(cr, incr_val, REF_CSY, "cy*", "y", pmodel); if ( val > 0.1 ) fail = 1; - incr_val = incr_frac * cx; - val = test_gradients(cr, incr_val, REF_CSX, "cx*", pmodel); + incr_val = incr_frac * az; + val = test_gradients(cr, incr_val, REF_ASZ, "az*", "z", pmodel); if ( val > 0.1 ) fail = 1; - incr_val = incr_frac * cy; - val = test_gradients(cr, incr_val, REF_CSY, "cy*", pmodel); + incr_val = incr_frac * bz; + val = test_gradients(cr, incr_val, REF_BSZ, "bz*", "z", pmodel); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * cz; - val = test_gradients(cr, incr_val, REF_CSZ, "cz*", pmodel); + val = test_gradients(cr, incr_val, REF_CSZ, "cz*", "z", pmodel); if ( val > 0.1 ) fail = 1; } -- cgit v1.2.3 From 803336b21468042f64c27ca262248755e5ea77b9 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 8 Mar 2013 10:53:58 +0100 Subject: Add script to plot gradients --- tests/plot_gradients | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) create mode 100755 tests/plot_gradients diff --git a/tests/plot_gradients b/tests/plot_gradients new file mode 100755 index 00000000..b5800d07 --- /dev/null +++ b/tests/plot_gradients @@ -0,0 +1,28 @@ +#!/bin/sh + +gnuplot -persist << EOF + +set key bottom right + +set xlabel "Calculated gradient" +set ylabel "Observed gradient" + + plot "gradient-test-x.dat" using 1:2 w p lc 1 pt 1 title "x" +replot "gradient-test-y.dat" using 1:2 w p lc 2 pt 1 title "y" +replot "gradient-test-z.dat" using 1:2 w p lc 3 pt 1 title "z" + +EOF + +gnuplot -persist << EOF +set key bottom right +set xlabel "Calculated gradient" +set ylabel "Observed gradient" +plot "gradient-test-R.dat" using 1:2 w p lc 1 pt 1 title "profile radius" +EOF + +gnuplot -persist << EOF +set key bottom right +set xlabel "Calculated gradient" +set ylabel "Observed gradient" +plot "gradient-test-div.dat" using 1:2 w p lc 1 pt 1 title "divergence" +EOF -- cgit v1.2.3 From 2d3209e0ad03ea6bc72a91f0bf34c7d5a3ba3732 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 8 Mar 2013 15:05:52 +0100 Subject: pr_gradient_check: Make sure list index keeps step with reflection list --- tests/pr_gradient_check.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/pr_gradient_check.c b/tests/pr_gradient_check.c index 0dd5466b..1ca3a424 100644 --- a/tests/pr_gradient_check.c +++ b/tests/pr_gradient_check.c @@ -295,6 +295,7 @@ static double test_gradients(Crystal *cr, double incr_val, int refine, if ( !valid[i] ) { n_invalid++; + i++; } else { double r1, r2, p; @@ -303,6 +304,7 @@ static double test_gradients(Crystal *cr, double incr_val, int refine, grad1 = (vals[1][i] - vals[0][i]) / incr_val; grad2 = (vals[2][i] - vals[1][i]) / incr_val; grad = (grad1 + grad2) / 2.0; + i++; cgrad = gradient(cr, refine, refl, pmodel); @@ -350,8 +352,6 @@ static double test_gradients(Crystal *cr, double incr_val, int refine, } - i++; - } STATUS("%3s: %3i within 5%%, %3i outside, %3i nan, %3i invalid, " -- cgit v1.2.3 From 698179b8bf305d0794d0bd205243c17682b7dc4a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 8 Mar 2013 15:23:15 +0100 Subject: pr_gradient_check: Improved checking --- tests/pr_gradient_check.c | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/tests/pr_gradient_check.c b/tests/pr_gradient_check.c index 1ca3a424..8896e58f 100644 --- a/tests/pr_gradient_check.c +++ b/tests/pr_gradient_check.c @@ -320,7 +320,7 @@ static double test_gradients(Crystal *cr, double incr_val, int refine, vec2[n_line] = grad; n_line++; - if ( fabs(cgrad) < 5e-8 ) { + if ( (fabs(cgrad) < 5e-8) && (fabs(grad) < 5e-8) ) { n_small++; continue; } @@ -328,14 +328,15 @@ static double test_gradients(Crystal *cr, double incr_val, int refine, total += fabs(cgrad - grad); ntot++; - if ( !within_tolerance(grad, cgrad, 10.0) ) + if ( !within_tolerance(grad, cgrad, 10.0) + || !within_tolerance(cgrad, grad, 10.0) ) { - //STATUS("!- %s %3i %3i %3i" - // " %10.2Le %10.2e ratio = %5.2Lf" - // " %10.2e %10.2e\n", - // str, h, k, l, grad, cgrad, cgrad/grad, - // r1, r2); + STATUS("!- %s %3i %3i %3i" + " %10.2Le %10.2e ratio = %5.2Lf" + " %10.2e %10.2e\n", + str, h, k, l, grad, cgrad, cgrad/grad, + r1, r2); n_bad++; } else { -- cgit v1.2.3 From c5ededeb46068dbb8f3d1bcb465611b4b4d78b04 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 8 Mar 2013 15:26:57 +0100 Subject: pr_gradient_check: Add --quiet, tidy up --- tests/pr_gradient_check.c | 66 +++++++++++++++++++++++++++++++---------------- 1 file changed, 44 insertions(+), 22 deletions(-) diff --git a/tests/pr_gradient_check.c b/tests/pr_gradient_check.c index 8896e58f..48719914 100644 --- a/tests/pr_gradient_check.c +++ b/tests/pr_gradient_check.c @@ -222,7 +222,7 @@ static void calc_either_side(Crystal *cr, double incr_val, static double test_gradients(Crystal *cr, double incr_val, int refine, const char *str, const char *file, - PartialityModel pmodel) + PartialityModel pmodel, int quiet) { Reflection *refl; RefListIterator *iter; @@ -328,15 +328,17 @@ static double test_gradients(Crystal *cr, double incr_val, int refine, total += fabs(cgrad - grad); ntot++; - if ( !within_tolerance(grad, cgrad, 10.0) - || !within_tolerance(cgrad, grad, 10.0) ) + if ( !within_tolerance(grad, cgrad, 5.0) + || !within_tolerance(cgrad, grad, 5.0) ) { - STATUS("!- %s %3i %3i %3i" - " %10.2Le %10.2e ratio = %5.2Lf" - " %10.2e %10.2e\n", - str, h, k, l, grad, cgrad, cgrad/grad, - r1, r2); + if ( !quiet ) { + STATUS("!- %s %3i %3i %3i" + " %10.2Le %10.2e ratio = %5.2Lf" + " %10.2e %10.2e\n", + str, h, k, l, grad, cgrad, + cgrad/grad, r1, r2); + } n_bad++; } else { @@ -379,6 +381,19 @@ int main(int argc, char *argv[]) int i; const PartialityModel pmodel = PMODEL_SPHERE; int fail = 0; + int quiet = 0; + + if ( argc == 2 ) { + if ( strcmp(argv[1], "--quiet") == 0 ) { + quiet = 1; + } else { + ERROR("Syntax: %s [--quiet]\n", argv[0]); + return 1; + } + } else if ( argc != 1 ) { + ERROR("Syntax: %s [--quiet]\n", argv[0]); + return 1; + } image.width = 1024; image.height = 1024; @@ -421,44 +436,51 @@ int main(int argc, char *argv[]) incr_val = incr_frac * image.div; val = test_gradients(cr, incr_val, REF_DIV, "div", "div", - pmodel); + pmodel, quiet); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * crystal_get_profile_radius(cr); - val = test_gradients(cr, incr_val, REF_R, "R", "R", pmodel); + val = test_gradients(cr, incr_val, REF_R, "R", "R", pmodel, + quiet); if ( val > 0.1 ) fail = 1; - //incr_val = incr_frac * image.profile_radius; - //val += test_gradients(&image, incr_val, REF_R, "rad"); - incr_val = incr_frac * ax; - val = test_gradients(cr, incr_val, REF_ASX, "ax*", "x", pmodel); + val = test_gradients(cr, incr_val, REF_ASX, "ax*", "x", pmodel, + quiet); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * bx; - val = test_gradients(cr, incr_val, REF_BSX, "bx*", "x", pmodel); + val = test_gradients(cr, incr_val, REF_BSX, "bx*", "x", pmodel, + quiet); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * cx; - val = test_gradients(cr, incr_val, REF_CSX, "cx*", "x", pmodel); + val = test_gradients(cr, incr_val, REF_CSX, "cx*", "x", pmodel, + quiet); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * ay; - val = test_gradients(cr, incr_val, REF_ASY, "ay*", "y", pmodel); + val = test_gradients(cr, incr_val, REF_ASY, "ay*", "y", pmodel, + quiet); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * by; - val = test_gradients(cr, incr_val, REF_BSY, "by*", "y", pmodel); + val = test_gradients(cr, incr_val, REF_BSY, "by*", "y", pmodel, + quiet); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * cy; - val = test_gradients(cr, incr_val, REF_CSY, "cy*", "y", pmodel); + val = test_gradients(cr, incr_val, REF_CSY, "cy*", "y", pmodel, + quiet); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * az; - val = test_gradients(cr, incr_val, REF_ASZ, "az*", "z", pmodel); + val = test_gradients(cr, incr_val, REF_ASZ, "az*", "z", pmodel, + quiet); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * bz; - val = test_gradients(cr, incr_val, REF_BSZ, "bz*", "z", pmodel); + val = test_gradients(cr, incr_val, REF_BSZ, "bz*", "z", pmodel, + quiet); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * cz; - val = test_gradients(cr, incr_val, REF_CSZ, "cz*", "z", pmodel); + val = test_gradients(cr, incr_val, REF_CSZ, "cz*", "z", pmodel, + quiet); if ( val > 0.1 ) fail = 1; } -- cgit v1.2.3 From d3b353ad71a4d53e108d8e1dfeaff09ebaf14e17 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 12 Mar 2013 10:52:14 +0100 Subject: Fix div gradients --- src/post-refinement.c | 56 +++++++++++++++++++-------------------------------- 1 file changed, 21 insertions(+), 35 deletions(-) diff --git a/src/post-refinement.c b/src/post-refinement.c index 9e63736f..32a93f38 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -91,8 +91,7 @@ static double partiality_rgradient(double r, double profile_radius) double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) { double ds, azi; - double nom, den; - double g; + double glow, ghigh; double asx, asy, asz; double bsx, bsy, bsz; double csx, csy, csz; @@ -142,12 +141,15 @@ double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) azi = atan2(yl, xl); /* Calculate the gradient of partiality wrt excitation error. */ - g = 0.0; if ( clamp_low == 0 ) { - g -= partiality_gradient(rlow, r); + glow = partiality_gradient(rlow, r); + } else { + glow = 0.0; } if ( clamp_high == 0 ) { - g += partiality_gradient(rhigh, r); + ghigh = partiality_gradient(rhigh, r); + } else { + ghigh = 0.0; } /* For many gradients, just multiply the above number by the gradient @@ -155,57 +157,41 @@ double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) switch ( k ) { case REF_DIV : - gr = 0.0; - if ( clamp_low == 0 ) { - nom = sqrt(2.0) * ds * sin(image->div/2.0); - den = sqrt(1.0 - cos(image->div/2.0)); - gr -= (nom/den) * g; - } - if ( clamp_high == 0 ) { - nom = sqrt(2.0) * ds * sin(image->div/2.0); - den = sqrt(1.0 - cos(image->div/2.0)); - gr += (nom/den) * g; - } - if ( isnan(gr) ) gr = 0.0; /* FIXME: This isn't true (?) */ - return gr / 4.0; /* FIXME: Shameless fudge factor */ + /* Small angle approximation */ + return (ds*glow + ds*ghigh) / 2.0; case REF_R : - g = 0.0; - if ( clamp_low == 0 ) { - g += partiality_rgradient(rlow, r); - } - if ( clamp_high == 0 ) { - g += partiality_rgradient(rhigh, r); - } - return g; + gr = partiality_rgradient(rlow, r); + gr += partiality_rgradient(rhigh, r); + return gr; /* Cell parameters and orientation */ case REF_ASX : - return hs * sin(phi) * cos(azi) * g; + return hs * sin(phi) * cos(azi) * (ghigh-glow); case REF_BSX : - return ks * sin(phi) * cos(azi) * g; + return ks * sin(phi) * cos(azi) * (ghigh-glow); case REF_CSX : - return ls * sin(phi) * cos(azi) * g; + return ls * sin(phi) * cos(azi) * (ghigh-glow); case REF_ASY : - return hs * sin(phi) * sin(azi) * g; + return hs * sin(phi) * sin(azi) * (ghigh-glow); case REF_BSY : - return ks * sin(phi) * sin(azi) * g; + return ks * sin(phi) * sin(azi) * (ghigh-glow); case REF_CSY : - return ls * sin(phi) * sin(azi) * g; + return ls * sin(phi) * sin(azi) * (ghigh-glow); case REF_ASZ : - return hs * cos(phi) * g; + return hs * cos(phi) * (ghigh-glow); case REF_BSZ : - return ks * cos(phi) * g; + return ks * cos(phi) * (ghigh-glow); case REF_CSZ : - return ls * cos(phi) * g; + return ls * cos(phi) * (ghigh-glow); } -- cgit v1.2.3 From 1acd2a02602a98a6286553868ed770427087c021 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 12 Mar 2013 12:26:14 +0100 Subject: Fix radius gradient --- src/post-refinement.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/post-refinement.c b/src/post-refinement.c index 32a93f38..9e4649a2 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -81,7 +81,7 @@ static double partiality_rgradient(double r, double profile_radius) dpdq = 6.0*(q-pow(q, 2.0)); /* dq/drad */ - dqdrad = 0.5 * (1.0 - r * pow(profile_radius, -2.0)); + dqdrad = -0.5 * r * pow(profile_radius, -2.0); return dpdq * dqdrad; } @@ -162,7 +162,7 @@ double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) case REF_R : gr = partiality_rgradient(rlow, r); - gr += partiality_rgradient(rhigh, r); + gr -= partiality_rgradient(rhigh, r); return gr; /* Cell parameters and orientation */ -- cgit v1.2.3 From 2b115b938d8bfba13f4193e329728b09a750b79f Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 13 Mar 2013 10:43:42 +0100 Subject: Add Lorentz factor --- doc/reference/libcrystfel/CrystFEL-sections.txt | 5 ++- libcrystfel/src/geometry.c | 19 ++++++++- libcrystfel/src/reflist.c | 51 ++++++++++++++++++++++--- libcrystfel/src/reflist.h | 9 +++-- src/partial_sim.c | 5 ++- 5 files changed, 76 insertions(+), 13 deletions(-) diff --git a/doc/reference/libcrystfel/CrystFEL-sections.txt b/doc/reference/libcrystfel/CrystFEL-sections.txt index ff8001cf..09267538 100644 --- a/doc/reference/libcrystfel/CrystFEL-sections.txt +++ b/doc/reference/libcrystfel/CrystFEL-sections.txt @@ -21,10 +21,11 @@ next_found_refl get_excitation_error get_detector_pos get_partiality +get_lorentz +get_partial get_indices get_symmetric_indices get_intensity -get_partial get_scalable get_refinable get_redundancy @@ -34,6 +35,8 @@ get_temp1 get_temp2 set_detector_pos +set_partiality +set_lorentz set_partial set_intensity set_scalable diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index 450ebfaf..7a553f82 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -134,6 +134,10 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, Reflection *refl; double cet, cez; double pr; + double L; + + /* Don't predict 000 */ + if ( abs(h)+abs(k)+abs(l) == 0 ) return NULL; pr = crystal_get_profile_radius(cryst); @@ -160,6 +164,19 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, && (fabs(rlow) > pr) && (fabs(rhigh) > pr) ) return NULL; + if ( rlow < rhigh ) { + ERROR("Reflection with rlow < rhigh!\n"); + ERROR("%3i %3i %3i rlow = %e, rhigh = %e\n", + h, k, l, rlow, rhigh); + ERROR("div = %e\n", image->div); + return NULL; + } + + /* Lorentz factor is determined direction from the r values, before + * clamping. The multiplication by the profile radius is to make the + * correction factor vaguely near 1. */ + L = pr / (rlow - rhigh); + /* If the "lower" Ewald sphere is a long way away, use the * position at which the Ewald sphere would just touch the * reflection. @@ -185,7 +202,6 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, rhigh = +pr; clamp_high = +1; } - assert(clamp_low >= clamp_high); /* Calculate partiality */ part = partiality(rlow, rhigh, pr); @@ -198,6 +214,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, refl = reflection_new(h, k, l); set_detector_pos(refl, 0.0, xda, yda); set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high); + set_lorentz(refl, L); set_symmetric_indices(refl, h, k, l); set_redundancy(refl, 1); diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c index 64154774..b3b9f85b 100644 --- a/libcrystfel/src/reflist.c +++ b/libcrystfel/src/reflist.c @@ -3,11 +3,11 @@ * * Fast reflection/peak list * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2011-2012 Thomas White + * 2011-2013 Thomas White * * This file is part of CrystFEL. * @@ -70,6 +70,7 @@ struct _refldata { double r1; /* First excitation error */ double r2; /* Second excitation error */ double p; /* Partiality */ + double L; /* Lorentz factor */ int clamp1; /* Clamp status for r1 */ int clamp2; /* Clamp status for r2 */ @@ -388,7 +389,7 @@ void get_symmetric_indices(const Reflection *refl, * get_partiality: * @refl: A %Reflection * - * Returns: The partiality of the reflection. + * Returns: The partiality of the reflection. See get_lorentz(). **/ double get_partiality(const Reflection *refl) { @@ -396,6 +397,19 @@ double get_partiality(const Reflection *refl) } +/** + * get_lorentz: + * @refl: A %Reflection + * + * Returns: The Lorentz factor for the reflection. To "scale up" a partial + * reflection, divide by this multiplied by the partiality. + **/ +double get_lorentz(const Reflection *refl) +{ + return refl->data.L; +} + + /** * get_intensity: * @refl: A %Reflection @@ -599,13 +613,38 @@ void set_partial(Reflection *refl, double r1, double r2, double p, } +/** + * set_intensity: + * @refl: A %Reflection + * @p: The partiality for the reflection. + * + * Set the partiality for the reflection. See set_lorentz(). + **/ +void set_partiality(Reflection *refl, double p) +{ + refl->data.p = p; +} + +/** + * set_lorentz: + * @refl: A %Reflection + * @L: The Lorentz factor for the reflection. + * + * Set the Lorentz factor for the reflection. To "scale up" a partial + * reflection, divide by this multiplied by the partiality. + **/ +void set_lorentz(Reflection *refl, double L) +{ + refl->data.L = L; +} + + /** * set_intensity: * @refl: A %Reflection * @intensity: The intensity for the reflection. * - * Set the intensity for the reflection. Note that retrieval is done with - * get_intensity(). + * Set the intensity for the reflection. **/ void set_intensity(Reflection *refl, double intensity) { diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h index 246ef885..7969235c 100644 --- a/libcrystfel/src/reflist.h +++ b/libcrystfel/src/reflist.h @@ -3,11 +3,11 @@ * * Fast reflection/peak list * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2011-2012 Thomas White + * 2011-2013 Thomas White * * This file is part of CrystFEL. * @@ -78,6 +78,7 @@ extern Reflection *next_found_refl(Reflection *refl); extern double get_excitation_error(const Reflection *refl); extern void get_detector_pos(const Reflection *refl, double *fs, double *ss); extern double get_partiality(const Reflection *refl); +extern double get_lorentz(const Reflection *refl); extern void get_indices(const Reflection *refl, signed int *h, signed int *k, signed int *l); extern void get_symmetric_indices(const Reflection *refl, @@ -100,6 +101,8 @@ extern void set_detector_pos(Reflection *refl, double exerr, double fs, double ss); extern void set_partial(Reflection *refl, double r1, double r2, double p, double clamp_low, double clamp_high); +extern void set_partiality(Reflection *refl, double p); +extern void set_lorentz(Reflection *refl, double L); extern void set_intensity(Reflection *refl, double intensity); extern void set_scalable(Reflection *refl, int scalable); extern void set_refinable(Reflection *refl, int refinable); diff --git a/src/partial_sim.c b/src/partial_sim.c index 7c93da1b..32e6666a 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -100,12 +100,13 @@ static void calculate_partials(Crystal *cr, { signed int h, k, l; Reflection *rfull; - double p, Ip, If; + double L, p, Ip, If; int bin; get_indices(refl, &h, &k, &l); get_asymm(sym, h, k, l, &h, &k, &l); p = get_partiality(refl); + L = get_lorentz(refl); pthread_mutex_lock(full_lock); rfull = find_refl(full, h, k, l); @@ -139,7 +140,7 @@ static void calculate_partials(Crystal *cr, } } - Ip = crystal_get_osf(cr) * p * If; + Ip = crystal_get_osf(cr) * L * p * If; res = resolution(crystal_get_cell(cr), h, k, l); bin = NBINS*2.0*res/max_q; -- cgit v1.2.3 From 060dcd17f88045995eb54e8f663e06f9144b13bf Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 13 Mar 2013 10:43:52 +0100 Subject: Fix comment --- libcrystfel/src/geometry.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h index 8b726ed7..2ac5982b 100644 --- a/libcrystfel/src/geometry.h +++ b/libcrystfel/src/geometry.h @@ -49,7 +49,7 @@ extern "C" { * space. * @PMODEL_UNITY : Set all all partialities and Lorentz factors to 1. * - * A %PartialityModelModel describes a geometrical model which can be used to + * A %PartialityModel describes a geometrical model which can be used to * calculate spot partialities and Lorentz correction factors. **/ typedef enum { -- cgit v1.2.3 From e4d9bbeedfc7664f4abfdee7ec7b565ae169e6d9 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 15 Mar 2013 11:02:41 +0100 Subject: tests/pr_gradient_check: Add --plot --- tests/pr_gradient_check.c | 72 +++++++++++++++++++++++++++++++---------------- 1 file changed, 47 insertions(+), 25 deletions(-) diff --git a/tests/pr_gradient_check.c b/tests/pr_gradient_check.c index 48719914..c7b70563 100644 --- a/tests/pr_gradient_check.c +++ b/tests/pr_gradient_check.c @@ -30,6 +30,7 @@ #include #include #include +#include #include #include @@ -222,7 +223,7 @@ static void calc_either_side(Crystal *cr, double incr_val, static double test_gradients(Crystal *cr, double incr_val, int refine, const char *str, const char *file, - PartialityModel pmodel, int quiet) + PartialityModel pmodel, int quiet, int plot) { Reflection *refl; RefListIterator *iter; @@ -269,8 +270,10 @@ static double test_gradients(Crystal *cr, double incr_val, int refine, calc_either_side(cr, incr_val, valid, vals, refine); - snprintf(tmp, 32, "gradient-test-%s.dat", file); - fh = fopen(tmp, "w"); + if ( plot ) { + snprintf(tmp, 32, "gradient-test-%s.dat", file); + fh = fopen(tmp, "w"); + } vec1 = malloc(nref*sizeof(double)); vec2 = malloc(nref*sizeof(double)); @@ -315,7 +318,10 @@ static double test_gradients(Crystal *cr, double incr_val, int refine, continue; } - fprintf(fh, "%e %Le\n", cgrad, grad); + if ( plot ) { + fprintf(fh, "%e %Le\n", cgrad, grad); + } + vec1[n_line] = cgrad; vec2[n_line] = grad; n_line++; @@ -359,7 +365,10 @@ static double test_gradients(Crystal *cr, double incr_val, int refine, STATUS("%3s: %3i within 5%%, %3i outside, %3i nan, %3i invalid, " "%3i small. ", str, n_good, n_bad, n_nan, n_invalid, n_small); - fclose(fh); + + if ( plot ) { + fclose(fh); + } cc = gsl_stats_correlation(vec1, 1, vec2, 1, n_line); STATUS("CC = %+f\n", cc); @@ -382,17 +391,30 @@ int main(int argc, char *argv[]) const PartialityModel pmodel = PMODEL_SPHERE; int fail = 0; int quiet = 0; + int plot = 0; + int c; + + const struct option longopts[] = { + {"quiet", 0, &quiet, 1}, + {"plot", 0, &plot, 1}, + {0, 0, NULL, 0} + }; + + while ((c = getopt_long(argc, argv, "", longopts, NULL)) != -1) { + switch (c) { + + case 0 : + break; + + case '?' : + break; + + default : + ERROR("Unhandled option '%c'\n", c); + break; - if ( argc == 2 ) { - if ( strcmp(argv[1], "--quiet") == 0 ) { - quiet = 1; - } else { - ERROR("Syntax: %s [--quiet]\n", argv[0]); - return 1; } - } else if ( argc != 1 ) { - ERROR("Syntax: %s [--quiet]\n", argv[0]); - return 1; + } image.width = 1024; @@ -436,51 +458,51 @@ int main(int argc, char *argv[]) incr_val = incr_frac * image.div; val = test_gradients(cr, incr_val, REF_DIV, "div", "div", - pmodel, quiet); + pmodel, quiet, plot); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * crystal_get_profile_radius(cr); val = test_gradients(cr, incr_val, REF_R, "R", "R", pmodel, - quiet); + quiet, plot); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * ax; val = test_gradients(cr, incr_val, REF_ASX, "ax*", "x", pmodel, - quiet); + quiet, plot); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * bx; val = test_gradients(cr, incr_val, REF_BSX, "bx*", "x", pmodel, - quiet); + quiet, plot); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * cx; val = test_gradients(cr, incr_val, REF_CSX, "cx*", "x", pmodel, - quiet); + quiet, plot); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * ay; val = test_gradients(cr, incr_val, REF_ASY, "ay*", "y", pmodel, - quiet); + quiet, plot); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * by; val = test_gradients(cr, incr_val, REF_BSY, "by*", "y", pmodel, - quiet); + quiet, plot); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * cy; val = test_gradients(cr, incr_val, REF_CSY, "cy*", "y", pmodel, - quiet); + quiet, plot); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * az; val = test_gradients(cr, incr_val, REF_ASZ, "az*", "z", pmodel, - quiet); + quiet, plot); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * bz; val = test_gradients(cr, incr_val, REF_BSZ, "bz*", "z", pmodel, - quiet); + quiet, plot); if ( val > 0.1 ) fail = 1; incr_val = incr_frac * cz; val = test_gradients(cr, incr_val, REF_CSZ, "cz*", "z", pmodel, - quiet); + quiet, plot); if ( val > 0.1 ) fail = 1; } -- cgit v1.2.3 From b31a273b8191bd0688d1382168ce9ba3500e6604 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 15 Mar 2013 11:04:17 +0100 Subject: tests/pr_gradient_check: Fix criterion --- tests/pr_gradient_check.c | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/tests/pr_gradient_check.c b/tests/pr_gradient_check.c index c7b70563..68a625a7 100644 --- a/tests/pr_gradient_check.c +++ b/tests/pr_gradient_check.c @@ -459,51 +459,51 @@ int main(int argc, char *argv[]) incr_val = incr_frac * image.div; val = test_gradients(cr, incr_val, REF_DIV, "div", "div", pmodel, quiet, plot); - if ( val > 0.1 ) fail = 1; + if ( val < 0.99 ) fail = 1; incr_val = incr_frac * crystal_get_profile_radius(cr); val = test_gradients(cr, incr_val, REF_R, "R", "R", pmodel, quiet, plot); - if ( val > 0.1 ) fail = 1; + if ( val < 0.99 ) fail = 1; incr_val = incr_frac * ax; val = test_gradients(cr, incr_val, REF_ASX, "ax*", "x", pmodel, quiet, plot); - if ( val > 0.1 ) fail = 1; + if ( val < 0.99 ) fail = 1; incr_val = incr_frac * bx; val = test_gradients(cr, incr_val, REF_BSX, "bx*", "x", pmodel, quiet, plot); - if ( val > 0.1 ) fail = 1; + if ( val < 0.99 ) fail = 1; incr_val = incr_frac * cx; val = test_gradients(cr, incr_val, REF_CSX, "cx*", "x", pmodel, quiet, plot); - if ( val > 0.1 ) fail = 1; + if ( val < 0.99 ) fail = 1; incr_val = incr_frac * ay; val = test_gradients(cr, incr_val, REF_ASY, "ay*", "y", pmodel, quiet, plot); - if ( val > 0.1 ) fail = 1; + if ( val < 0.99 ) fail = 1; incr_val = incr_frac * by; val = test_gradients(cr, incr_val, REF_BSY, "by*", "y", pmodel, quiet, plot); - if ( val > 0.1 ) fail = 1; + if ( val < 0.99 ) fail = 1; incr_val = incr_frac * cy; val = test_gradients(cr, incr_val, REF_CSY, "cy*", "y", pmodel, quiet, plot); - if ( val > 0.1 ) fail = 1; + if ( val < 0.99 ) fail = 1; incr_val = incr_frac * az; val = test_gradients(cr, incr_val, REF_ASZ, "az*", "z", pmodel, quiet, plot); - if ( val > 0.1 ) fail = 1; + if ( val < 0.99 ) fail = 1; incr_val = incr_frac * bz; val = test_gradients(cr, incr_val, REF_BSZ, "bz*", "z", pmodel, quiet, plot); - if ( val > 0.1 ) fail = 1; + if ( val < 0.99 ) fail = 1; incr_val = incr_frac * cz; val = test_gradients(cr, incr_val, REF_CSZ, "cz*", "z", pmodel, quiet, plot); - if ( val > 0.1 ) fail = 1; + if ( val < 0.99 ) fail = 1; } -- cgit v1.2.3 From 454b01b732b98a0b1cddbd85f91df64560416494 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 15 Mar 2013 11:20:31 +0100 Subject: Use Lorentz factors when merging reflections --- libcrystfel/src/geometry.c | 4 +++- src/hrs-scaling.c | 25 +++++++++++++++---------- 2 files changed, 18 insertions(+), 11 deletions(-) diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index 7a553f82..f457fc4f 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -384,7 +384,7 @@ void update_partialities(Crystal *cryst, PartialityModel pmodel) refl = next_refl(refl, iter) ) { Reflection *vals; - double r1, r2, p, x, y; + double r1, r2, L, p, x, y; double xl, yl, zl; signed int h, k, l; int clamp1, clamp2; @@ -407,6 +407,8 @@ void update_partialities(Crystal *cryst, PartialityModel pmodel) /* Transfer partiality stuff */ get_partial(vals, &r1, &r2, &p, &clamp1, &clamp2); set_partial(refl, r1, r2, p, clamp1, clamp2); + L = get_lorentz(vals); + set_lorentz(refl, L); /* Transfer detector location */ get_detector_pos(vals, &x, &y); diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index 0498a532..f60dfe25 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -3,11 +3,11 @@ * * Intensity scaling using generalised HRS target function * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2010-2012 Thomas White + * 2010-2013 Thomas White * * This file is part of CrystFEL. * @@ -109,6 +109,7 @@ static void run_scale_job(void *vwargs, int cookie) signed int h, k, l; double Ih, Ihl, esd; Reflection *r; + double L, p; if ( !get_scalable(refl) ) continue; @@ -125,9 +126,11 @@ static void run_scale_job(void *vwargs, int cookie) Ih = get_intensity(r); } - /* FIXME: This isn't correct */ - Ihl = get_intensity(refl) / get_partiality(refl); - esd = get_esd_intensity(refl) / get_partiality(refl); + /* If you change this, be sure to change run_merge_job() too */ + p = get_partiality(refl); + L = get_lorentz(refl); + Ihl = get_intensity(refl) / (L*p); + esd = get_esd_intensity(refl) / (L*p); num += Ih * (Ihl/osf) / pow(esd/osf, 2.0); den += pow(Ih, 2.0)/pow(esd/osf, 2.0); @@ -228,7 +231,7 @@ static void run_merge_job(void *vwargs, int cookie) signed int h, k, l; double num, den; int red; - double Ihl, esd, pcalc; + double Ihl, esd, L, p; if ( !get_scalable(refl) ) continue; @@ -251,9 +254,11 @@ static void run_merge_job(void *vwargs, int cookie) red = get_redundancy(f); } - pcalc = get_partiality(refl); - Ihl = get_intensity(refl) / pcalc; - esd = get_esd_intensity(refl) / pcalc; + /* If you change this, be sure to change run_scale_job() too */ + p = get_partiality(refl); + L = get_lorentz(refl); + Ihl = get_intensity(refl) / (L*p); + esd = get_esd_intensity(refl) / (L*p); num += (Ihl/G) / pow(esd/G, 2.0); den += 1.0 / pow(esd/G, 2.0); -- cgit v1.2.3 From 101aa5f9f8290b1d180d04d6ab5ff7eac934c148 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 15 Mar 2013 11:30:47 +0100 Subject: partial_sim: Add --osf-stddev --- doc/man/partial_sim.1 | 6 ++++++ src/partial_sim.c | 34 +++++++++++++++++++++++++++++----- 2 files changed, 35 insertions(+), 5 deletions(-) diff --git a/doc/man/partial_sim.1 b/doc/man/partial_sim.1 index 5b417978..9f3436ea 100644 --- a/doc/man/partial_sim.1 +++ b/doc/man/partial_sim.1 @@ -77,6 +77,12 @@ Add random values with a flat distribution to the components of the reciprocal l .PD Save a table of values to \fIfilename\fR containing, in resolution shells, the following columns: resolution shell centre in nm^-1, number of reflections in shell, mean partiality, maximum partiality. +.PD 0 +.B +.IP "\fB--osf-stddev=\fR\fIval\fR +.PD +Set the standard deviation of the distribution of overall scaling factors to \fIval\fR. The distribution will be cut at zero, i.e. negative or zero scaling factors are not allowed. The distribution will be centered on 1. The default is \fB--osf-stddev=2.0\fR. + .SH AUTHOR This page was written by Thomas White. diff --git a/src/partial_sim.c b/src/partial_sim.c index 32e6666a..615c66a3 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -3,11 +3,11 @@ * * Generate partials for testing scaling * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2011-2012 Thomas White + * 2011-2013 Thomas White * * This file is part of CrystFEL. * @@ -184,6 +184,7 @@ static void show_help(const char *s) " -c, --cnoise= Add random noise, with a flat distribution, to the\n" " reciprocal lattice vector components given in the\n" " stream, with maximum error +/- percent.\n" +" --osf-stddev= Set the standard deviation of the scaling factors.\n" "\n" ); } @@ -202,6 +203,7 @@ struct queue_args int random_intensities; UnitCell *cell; double cnoise; + double osf_stddev; struct image *template_image; double max_q; @@ -255,6 +257,7 @@ static void run_job(void *vwargs, int cookie) int i; Crystal *cr; RefList *reflections; + double osf; cr = crystal_new(); if ( cr == NULL ) { @@ -264,7 +267,10 @@ static void run_job(void *vwargs, int cookie) wargs->crystal = cr; crystal_set_image(cr, &wargs->image); - crystal_set_osf(cr, gaussian_noise(1.0, 0.3)); + do { + osf = gaussian_noise(1.0, qargs->osf_stddev); + } while ( osf <= 0.0 ); + crystal_set_osf(cr, osf); crystal_set_profile_radius(cr, wargs->image.beam->profile_radius); /* Set up a random orientation */ @@ -344,6 +350,7 @@ int main(int argc, char *argv[]) int i; FILE *fh; char *phist_file = NULL; + double osf_stddev = 2.0; /* Long options */ const struct option longopts[] = { @@ -355,8 +362,11 @@ int main(int argc, char *argv[]) {"geometry", 1, NULL, 'g'}, {"symmetry", 1, NULL, 'y'}, {"save-random", 1, NULL, 'r'}, - {"pgraph", 1, NULL, 2}, {"cnoise", 1, NULL, 'c'}, + + {"pgraph", 1, NULL, 2}, + {"osf-stddev", 1, NULL, 3}, + {0, 0, NULL, 0} }; @@ -418,6 +428,19 @@ int main(int argc, char *argv[]) phist_file = strdup(optarg); break; + case 3 : + osf_stddev = strtod(optarg, &rval); + if ( *rval != '\0' ) { + ERROR("Invalid OSF standard deviation.\n"); + return 1; + } + if ( osf_stddev <= 0.0 ) { + ERROR("Invalid OSF standard deviation."); + ERROR(" (must be positive).\n"); + return 1; + } + break; + case 0 : break; @@ -549,6 +572,7 @@ int main(int argc, char *argv[]) qargs.template_image = ℑ qargs.stream = stream; qargs.cnoise = cnoise; + qargs.osf_stddev = osf_stddev; qargs.max_q = largest_q(&image); for ( i=0; i Date: Fri, 15 Mar 2013 11:44:22 +0100 Subject: partial_sim: Add --full-stddev and --noise-stddev --- doc/man/partial_sim.1 | 18 +++++++++++++++--- src/partial_sim.c | 49 ++++++++++++++++++++++++++++++++++++++++++++----- 2 files changed, 59 insertions(+), 8 deletions(-) diff --git a/doc/man/partial_sim.1 b/doc/man/partial_sim.1 index 9f3436ea..7321b9ff 100644 --- a/doc/man/partial_sim.1 +++ b/doc/man/partial_sim.1 @@ -1,7 +1,7 @@ .\" .\" partial_sim man page .\" -.\" Copyright © 2012 Thomas White +.\" Copyright © 2012-2013 Thomas White .\" .\" Part of CrystFEL - crystallography with a FEL .\" @@ -81,7 +81,19 @@ Save a table of values to \fIfilename\fR containing, in resolution shells, the f .B .IP "\fB--osf-stddev=\fR\fIval\fR .PD -Set the standard deviation of the distribution of overall scaling factors to \fIval\fR. The distribution will be cut at zero, i.e. negative or zero scaling factors are not allowed. The distribution will be centered on 1. The default is \fB--osf-stddev=2.0\fR. +Set the standard deviation of the distribution of overall scaling factors to \fIval\fR. The distribution will be cut at zero, i.e. negative or zero scaling factors are not allowed. The distribution will be Gaussian centered on 1. The default is \fB--osf-stddev=2.0\fR. + +.PD 0 +.B +.IP "\fB--full-stddev=\fR\fIval\fR +.PD +Set the standard deviation of the distribution of randomly generated full intensities to \fIval\fR. The distribution will be Gaussian, centered on zero, and the absolute value will be taken (i.e. there will be no negative values). The default is \fB--full-stddev=1000.0\fR. This option has no effect if you also use \fB-i\fR or \fR--input\fB. + +.PD 0 +.B +.IP "\fB--noise-stddev=\fR\fIval\fR +.PD +Set the standard deviation of the noise added to the partial intensities to \fIval\fR. The noise will be Gaussian, and the same for all reflections. The default is \fB--noise-stddev=20.0\fR. .SH AUTHOR This page was written by Thomas White. @@ -90,7 +102,7 @@ This page was written by Thomas White. Report bugs to , or visit . .SH COPYRIGHT AND DISCLAIMER -Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, a research centre of the Helmholtz Association. +Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, a research centre of the Helmholtz Association. .P partial_sim, and this manaul, are part of CrystFEL. .P diff --git a/src/partial_sim.c b/src/partial_sim.c index 615c66a3..4fcbc823 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -88,7 +88,8 @@ static void calculate_partials(Crystal *cr, int random_intensities, pthread_mutex_t *full_lock, unsigned long int *n_ref, double *p_hist, - double *p_max, double max_q) + double *p_max, double max_q, double full_stddev, + double noise_stddev) { Reflection *refl; RefListIterator *iter; @@ -121,7 +122,7 @@ static void calculate_partials(Crystal *cr, * thing under lock. */ pthread_mutex_lock(full_lock); rfull = add_refl(full, h, k, l); - If = fabs(gaussian_noise(0.0, 1000.0)); + If = fabs(gaussian_noise(0.0, full_stddev)); set_intensity(rfull, If); set_redundancy(rfull, 1); pthread_mutex_unlock(full_lock); @@ -153,10 +154,10 @@ static void calculate_partials(Crystal *cr, res, bin, p); } - Ip = gaussian_noise(Ip, 100.0); + Ip = gaussian_noise(Ip, noise_stddev); set_intensity(refl, Ip); - set_esd_intensity(refl, 100.0); + set_esd_intensity(refl, noise_stddev); } } @@ -185,6 +186,9 @@ static void show_help(const char *s) " reciprocal lattice vector components given in the\n" " stream, with maximum error +/- percent.\n" " --osf-stddev= Set the standard deviation of the scaling factors.\n" +" --full-stddev= Set the standard deviation of the randomly\n" +" generated full intensities, if not using -i.\n" +" --noise-stddev= Set the standard deviation of the noise.\n" "\n" ); } @@ -204,6 +208,8 @@ struct queue_args UnitCell *cell; double cnoise; double osf_stddev; + double full_stddev; + double noise_stddev; struct image *template_image; double max_q; @@ -291,7 +297,8 @@ static void run_job(void *vwargs, int cookie) qargs->sym, qargs->random_intensities, &qargs->full_lock, wargs->n_ref, wargs->p_hist, wargs->p_max, - qargs->max_q); + qargs->max_q, qargs->full_stddev, + qargs->noise_stddev); /* Give a slightly incorrect cell in the stream */ mess_up_cell(cr, qargs->cnoise); @@ -351,6 +358,8 @@ int main(int argc, char *argv[]) FILE *fh; char *phist_file = NULL; double osf_stddev = 2.0; + double full_stddev = 1000.0; + double noise_stddev = 20.0; /* Long options */ const struct option longopts[] = { @@ -366,6 +375,8 @@ int main(int argc, char *argv[]) {"pgraph", 1, NULL, 2}, {"osf-stddev", 1, NULL, 3}, + {"full-stddev", 1, NULL, 4}, + {"noise-stddev", 1, NULL, 5}, {0, 0, NULL, 0} }; @@ -441,6 +452,32 @@ int main(int argc, char *argv[]) } break; + case 4 : + full_stddev = strtod(optarg, &rval); + if ( *rval != '\0' ) { + ERROR("Invalid full standard deviation.\n"); + return 1; + } + if ( full_stddev <= 0.0 ) { + ERROR("Invalid full standard deviation."); + ERROR(" (must be positive).\n"); + return 1; + } + break; + + case 5 : + noise_stddev = strtod(optarg, &rval); + if ( *rval != '\0' ) { + ERROR("Invalid noise standard deviation.\n"); + return 1; + } + if ( noise_stddev <= 0.0 ) { + ERROR("Invalid noise standard deviation."); + ERROR(" (must be positive).\n"); + return 1; + } + break; + case 0 : break; @@ -573,6 +610,8 @@ int main(int argc, char *argv[]) qargs.stream = stream; qargs.cnoise = cnoise; qargs.osf_stddev = osf_stddev; + qargs.full_stddev = full_stddev; + qargs.noise_stddev = noise_stddev; qargs.max_q = largest_q(&image); for ( i=0; i Date: Fri, 15 Mar 2013 11:46:40 +0100 Subject: partial_sim: Standard deviations can be zero --- src/partial_sim.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/partial_sim.c b/src/partial_sim.c index 4fcbc823..87d30f39 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -445,7 +445,7 @@ int main(int argc, char *argv[]) ERROR("Invalid OSF standard deviation.\n"); return 1; } - if ( osf_stddev <= 0.0 ) { + if ( osf_stddev < 0.0 ) { ERROR("Invalid OSF standard deviation."); ERROR(" (must be positive).\n"); return 1; @@ -458,7 +458,7 @@ int main(int argc, char *argv[]) ERROR("Invalid full standard deviation.\n"); return 1; } - if ( full_stddev <= 0.0 ) { + if ( full_stddev < 0.0 ) { ERROR("Invalid full standard deviation."); ERROR(" (must be positive).\n"); return 1; @@ -471,7 +471,7 @@ int main(int argc, char *argv[]) ERROR("Invalid noise standard deviation.\n"); return 1; } - if ( noise_stddev <= 0.0 ) { + if ( noise_stddev < 0.0 ) { ERROR("Invalid noise standard deviation."); ERROR(" (must be positive).\n"); return 1; -- cgit v1.2.3 From cef0a71eb385773fa2290dfc99de225948fef06a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 15 Mar 2013 15:53:39 +0100 Subject: partialator: Add --model= --- libcrystfel/src/geometry.c | 23 +++++++++++ src/hrs-scaling.c | 101 ++++++++++++++++++++++++++++++++++----------- src/hrs-scaling.h | 9 ++-- src/partialator.c | 17 +++++--- 4 files changed, 116 insertions(+), 34 deletions(-) diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index f457fc4f..cdc62ab5 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -362,6 +362,21 @@ RefList *select_intersections(struct image *image, Crystal *cryst) } +static void set_unity_partialities(Crystal *cryst) +{ + Reflection *refl; + RefListIterator *iter; + + for ( refl = first_refl(crystal_get_reflections(cryst), &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + set_partiality(refl, 1.0); + set_lorentz(refl, 1.0); + } +} + + /* Calculate partialities and apply them to the image's reflections */ void update_partialities(Crystal *cryst, PartialityModel pmodel) { @@ -373,6 +388,14 @@ void update_partialities(Crystal *cryst, PartialityModel pmodel) double csx, csy, csz; struct image *image = crystal_get_image(cryst); + if ( pmodel == PMODEL_UNITY ) { + /* It isn't strictly necessary to set the partialities to 1, + * because the scaling stuff will just a correction factor of + * 1 anyway. This is just to help things make sense. */ + set_unity_partialities(cryst); + return; + } + cell_get_reciprocal(crystal_get_cell(cryst), &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index f60dfe25..3a51e6d9 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -60,6 +60,7 @@ struct scale_queue_args Crystal **crystals; int n_started; double max_shift; + PartialityModel pmodel; }; @@ -68,6 +69,7 @@ struct scale_worker_args Crystal *crystal; double shift; RefList *reference; + PartialityModel pmodel; }; @@ -78,6 +80,7 @@ static void *create_scale_job(void *vqargs) wargs = malloc(sizeof(struct scale_worker_args)); wargs->reference = qargs->reference; + wargs->pmodel = qargs->pmodel; wargs->crystal = qargs->crystals[qargs->n_started++]; @@ -109,7 +112,7 @@ static void run_scale_job(void *vwargs, int cookie) signed int h, k, l; double Ih, Ihl, esd; Reflection *r; - double L, p; + double corr; if ( !get_scalable(refl) ) continue; @@ -126,11 +129,21 @@ static void run_scale_job(void *vwargs, int cookie) Ih = get_intensity(r); } - /* If you change this, be sure to change run_merge_job() too */ - p = get_partiality(refl); - L = get_lorentz(refl); - Ihl = get_intensity(refl) / (L*p); - esd = get_esd_intensity(refl) / (L*p); + /* If you change this, be sure to also change + * run_merge_job() and run_esd_job(). */ + switch ( wargs->pmodel ) { + + case PMODEL_UNITY : + corr = 1.0; + break; + + case PMODEL_SPHERE : + corr = get_partiality(refl) * get_lorentz(refl); + break; + + } + Ihl = get_intensity(refl) / corr; + esd = get_esd_intensity(refl) / corr; num += Ih * (Ihl/osf) / pow(esd/osf, 2.0); den += pow(Ih, 2.0)/pow(esd/osf, 2.0); @@ -160,7 +173,7 @@ static void finalise_scale_job(void *vqargs, void *vwargs) static double iterate_scale(Crystal **crystals, int n, RefList *reference, - int n_threads) + int n_threads, PartialityModel pmodel) { struct scale_queue_args qargs; @@ -170,6 +183,7 @@ static double iterate_scale(Crystal **crystals, int n, RefList *reference, qargs.n_started = 0; qargs.crystals = crystals; qargs.max_shift = 0.0; + qargs.pmodel = pmodel; run_threads(n_threads, run_scale_job, create_scale_job, finalise_scale_job, &qargs, n, 0, 0, 0); @@ -184,6 +198,7 @@ struct merge_queue_args pthread_mutex_t full_lock; Crystal **crystals; int n_started; + PartialityModel pmodel; }; @@ -192,6 +207,7 @@ struct merge_worker_args Crystal *crystal; RefList *full; pthread_mutex_t *full_lock; + PartialityModel pmodel; }; @@ -203,6 +219,7 @@ static void *create_merge_job(void *vqargs) wargs = malloc(sizeof(struct merge_worker_args)); wargs->full = qargs->full; wargs->full_lock = &qargs->full_lock; + wargs->pmodel = qargs->pmodel; wargs->crystal = qargs->crystals[qargs->n_started++]; @@ -231,7 +248,7 @@ static void run_merge_job(void *vwargs, int cookie) signed int h, k, l; double num, den; int red; - double Ihl, esd, L, p; + double Ihl, esd, corr; if ( !get_scalable(refl) ) continue; @@ -254,11 +271,22 @@ static void run_merge_job(void *vwargs, int cookie) red = get_redundancy(f); } - /* If you change this, be sure to change run_scale_job() too */ - p = get_partiality(refl); - L = get_lorentz(refl); - Ihl = get_intensity(refl) / (L*p); - esd = get_esd_intensity(refl) / (L*p); + /* If you change this, be sure to also change + * run_scale_job() and run_esd_job(). */ + switch ( wargs->pmodel ) { + + case PMODEL_UNITY : + corr = 1.0; + break; + + case PMODEL_SPHERE : + corr = get_partiality(refl) * get_lorentz(refl); + break;; + + } + + Ihl = get_intensity(refl) / corr; + esd = get_esd_intensity(refl) / corr; num += (Ihl/G) / pow(esd/G, 2.0); den += 1.0 / pow(esd/G, 2.0); @@ -278,7 +306,8 @@ static void finalise_merge_job(void *vqargs, void *vwargs) } -static RefList *lsq_intensities(Crystal **crystals, int n, int n_threads) +static RefList *lsq_intensities(Crystal **crystals, int n, int n_threads, + PartialityModel pmodel) { RefList *full; struct merge_queue_args qargs; @@ -290,6 +319,7 @@ static RefList *lsq_intensities(Crystal **crystals, int n, int n_threads) qargs.full = full; qargs.n_started = 0; qargs.crystals = crystals; + qargs.pmodel = pmodel; pthread_mutex_init(&qargs.full_lock, NULL); run_threads(n_threads, run_merge_job, create_merge_job, @@ -316,6 +346,7 @@ struct esd_queue_args RefList *full; Crystal **crystals; int n_started; + PartialityModel pmodel; }; @@ -323,6 +354,7 @@ struct esd_worker_args { Crystal *crystal; RefList *full; + PartialityModel pmodel; }; @@ -333,6 +365,7 @@ static void *create_esd_job(void *vqargs) wargs = malloc(sizeof(struct esd_worker_args)); wargs->full = qargs->full; + wargs->pmodel = qargs->pmodel; wargs->crystal = qargs->crystals[qargs->n_started++]; @@ -360,7 +393,7 @@ static void run_esd_job(void *vwargs, int cookie) Reflection *f; signed int h, k, l; double num; - double Ihl, Ih; + double Ihl, Ih, corr; if ( !get_scalable(refl) ) continue; @@ -372,8 +405,23 @@ static void run_esd_job(void *vwargs, int cookie) num = get_temp1(f); + /* If you change this, be sure to also change + * run_scale_job() and run_merge_job(). */ + switch ( wargs->pmodel ) { + + case PMODEL_UNITY : + corr = 1.0; + break; + + case PMODEL_SPHERE : + corr = get_partiality(refl) * get_lorentz(refl); + break;; + + } + + Ih = get_intensity(f); - Ihl = get_intensity(refl) / (get_partiality(refl) * G); + Ihl = get_intensity(refl) / (corr * G); num += pow(Ihl - Ih, 2.0); @@ -390,7 +438,7 @@ static void finalise_esd_job(void *vqargs, void *vwargs) static void calculate_esds(Crystal **crystals, int n, RefList *full, - int n_threads, int min_red) + int n_threads, int min_red, PartialityModel pmodel) { struct esd_queue_args qargs; Reflection *refl; @@ -399,6 +447,7 @@ static void calculate_esds(Crystal **crystals, int n, RefList *full, qargs.full = full; qargs.n_started = 0; qargs.crystals = crystals; + qargs.pmodel = pmodel; for ( refl = first_refl(full, &iter); refl != NULL; @@ -430,7 +479,7 @@ static void calculate_esds(Crystal **crystals, int n, RefList *full, /* Scale the stack of images */ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref, - int n_threads, int noscale) + int n_threads, int noscale, PartialityModel pmodel) { int i; double max_corr; @@ -440,14 +489,15 @@ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref, for ( i=0; i create an initial list to refine against */ if ( gref == NULL ) { - full = lsq_intensities(crystals, n, n_threads); + full = lsq_intensities(crystals, n, n_threads, pmodel); } /* Iterate */ @@ -463,14 +513,15 @@ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref, reference = full; } - max_corr = iterate_scale(crystals, n, reference, n_threads); + max_corr = iterate_scale(crystals, n, reference, n_threads, + pmodel); //STATUS("Scaling iteration %2i: max correction = %5.2f\n", // i+1, max_corr); /* No reference -> generate list for next iteration */ if ( gref == NULL ) { reflist_free(full); - full = lsq_intensities(crystals, n, n_threads); + full = lsq_intensities(crystals, n, n_threads, pmodel); } //show_scale_factors(images, n); @@ -480,10 +531,10 @@ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref, } while ( (max_corr > 0.01) && (i < MAX_CYCLES) ); if ( gref != NULL ) { - full = lsq_intensities(crystals, n, n_threads); + full = lsq_intensities(crystals, n, n_threads, pmodel); } /* else we already did it */ - calculate_esds(crystals, n, full, n_threads, min_redundancy); + calculate_esds(crystals, n, full, n_threads, min_redundancy, pmodel); return full; } diff --git a/src/hrs-scaling.h b/src/hrs-scaling.h index 80940347..8ae12380 100644 --- a/src/hrs-scaling.h +++ b/src/hrs-scaling.h @@ -3,11 +3,11 @@ * * Intensity scaling using generalised HRS target function * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2010-2012 Thomas White + * 2010-2013 Thomas White * * This file is part of CrystFEL. * @@ -37,10 +37,11 @@ #include "crystal.h" #include "reflist.h" +#include "geometry.h" extern RefList *scale_intensities(Crystal **crystals, int n, RefList *reference, int n_threads, - int noscale); + int noscale, PartialityModel pmodel); #endif /* HRS_SCALING_H */ diff --git a/src/partialator.c b/src/partialator.c index 079e078e..bdacd62b 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -79,6 +79,7 @@ static void show_help(const char *s) " -r, --reference= Refine images against reflections in ,\n" " instead of taking the mean of the intensity\n" " estimates.\n" +" -m, --model= Specify partiality model.\n" "\n" " -j Run analyses in parallel.\n"); } @@ -140,13 +141,18 @@ static void done_image(void *vqargs, void *task) static void refine_all(Crystal **crystals, int n_crystals, struct detector *det, - RefList *full, int nthreads) + RefList *full, int nthreads, PartialityModel pmodel) { struct refine_args task_defaults; struct queue_args qargs; + /* If the partiality model is "p=1", this refinement is really, really + * easy... */ + if ( pmodel == PMODEL_UNITY ) return; + task_defaults.full = full; task_defaults.crystal = NULL; + task_defaults.pmodel = pmodel; qargs.task_defaults = task_defaults; qargs.n_started = 0; @@ -329,7 +335,8 @@ int main(int argc, char *argv[]) {"iterations", 1, NULL, 'n'}, {"no-scale", 0, &noscale, 1}, {"reference", 1, NULL, 'r'}, - {"partiality", 1, NULL, 'm'}, + {"model", 1, NULL, 'm'}, + {0, 0, NULL, 0} }; @@ -574,7 +581,7 @@ int main(int argc, char *argv[]) STATUS("\nPerforming initial scaling.\n"); if ( noscale ) STATUS("Scale factors fixed at 1.\n"); full = scale_intensities(crystals, n_crystals, reference, - nthreads, noscale); + nthreads, noscale, pmodel); sr = sr_titlepage(crystals, n_crystals, "scaling-report.pdf", infile, cmdline); @@ -597,7 +604,7 @@ int main(int argc, char *argv[]) /* Refine the geometry of all patterns to get the best fit */ select_reflections_for_refinement(crystals, n_crystals, comp, have_reference); - refine_all(crystals, n_crystals, det, comp, nthreads); + refine_all(crystals, n_crystals, det, comp, nthreads, pmodel); nobs = 0; for ( j=0; j Date: Tue, 19 Mar 2013 14:33:45 +0100 Subject: Typo --- src/hrs-scaling.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index 3a51e6d9..fb1466cb 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -281,7 +281,7 @@ static void run_merge_job(void *vwargs, int cookie) case PMODEL_SPHERE : corr = get_partiality(refl) * get_lorentz(refl); - break;; + break; } -- cgit v1.2.3 From 3baecd0039864ac1e95fbc0b397fd203d02f4188 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 21 Mar 2013 16:24:04 +0100 Subject: s/osf/G/ --- src/hrs-scaling.c | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index fb1466cb..d91a36d9 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -97,8 +97,8 @@ static void run_scale_job(void *vwargs, int cookie) RefListIterator *iter; double num = 0.0; double den = 0.0; - double corr; - const double osf = crystal_get_osf(cr); + double g; + const double G = crystal_get_osf(cr); if ( crystal_get_user_flag(cr) ) { wargs->shift = 0.0; @@ -145,20 +145,19 @@ static void run_scale_job(void *vwargs, int cookie) Ihl = get_intensity(refl) / corr; esd = get_esd_intensity(refl) / corr; - num += Ih * (Ihl/osf) / pow(esd/osf, 2.0); - den += pow(Ih, 2.0)/pow(esd/osf, 2.0); + num += Ih * (Ihl/G) / pow(esd/G, 2.0); + den += pow(Ih, 2.0)/pow(esd/G, 2.0); } //num += image->osf / pow(SCALING_RESTRAINT, 2.0); //den += pow(image->osf, 2.0)/pow(SCALING_RESTRAINT, 2.0); - corr = num / den; - if ( !isnan(corr) && !isinf(corr) ) { - crystal_set_osf(cr, osf*corr); + g = num / den; + if ( !isnan(g) && !isinf(g) ) { + crystal_set_osf(cr, g*G); } - wargs->shift = fabs(corr-1.0); - + wargs->shift = fabs(g-1.0); } -- cgit v1.2.3 From 74970d1f8fc6476d6f5719beaa26cfbc843f8ada Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 21 Mar 2013 16:24:14 +0100 Subject: Spot bad partiality models (and shut the compiler up) --- src/hrs-scaling.c | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index d91a36d9..5acb247f 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -141,7 +141,13 @@ static void run_scale_job(void *vwargs, int cookie) corr = get_partiality(refl) * get_lorentz(refl); break; + default : + ERROR("Unrecognised partiality model!\n"); + abort(); + break; + } + Ihl = get_intensity(refl) / corr; esd = get_esd_intensity(refl) / corr; @@ -282,6 +288,11 @@ static void run_merge_job(void *vwargs, int cookie) corr = get_partiality(refl) * get_lorentz(refl); break; + default : + ERROR("Unrecognised partiality model!\n"); + abort(); + break; + } Ihl = get_intensity(refl) / corr; @@ -416,8 +427,12 @@ static void run_esd_job(void *vwargs, int cookie) corr = get_partiality(refl) * get_lorentz(refl); break;; - } + default : + ERROR("Unrecognised partiality model!\n"); + abort(); + break; + } Ih = get_intensity(f); Ihl = get_intensity(refl) / (corr * G); -- cgit v1.2.3 From 79ff2b85e3ef35fefe088bffbf50bfd195e7363d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 21 Mar 2013 16:24:39 +0100 Subject: Show number of scalable observations --- src/partialator.c | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/partialator.c b/src/partialator.c index bdacd62b..58025b96 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -554,6 +554,7 @@ int main(int argc, char *argv[]) display_progress(n_images, n_crystals); } while ( 1 ); + fprintf(stderr, "\n"); close_stream(st); @@ -576,9 +577,10 @@ int main(int argc, char *argv[]) } } + STATUS("%i scalable observations.\n", nobs); /* Make initial estimates */ - STATUS("\nPerforming initial scaling.\n"); + STATUS("Performing initial scaling.\n"); if ( noscale ) STATUS("Scale factors fixed at 1.\n"); full = scale_intensities(crystals, n_crystals, reference, nthreads, noscale, pmodel); -- cgit v1.2.3 From bdf3a0acbe270ea081bbad290f55767991c71c9c Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 21 Mar 2013 16:26:35 +0100 Subject: Add tests/partialator_merge_check_1 --- tests/partialator_merge_check_1 | 76 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 76 insertions(+) create mode 100755 tests/partialator_merge_check_1 diff --git a/tests/partialator_merge_check_1 b/tests/partialator_merge_check_1 new file mode 100755 index 00000000..633cb815 --- /dev/null +++ b/tests/partialator_merge_check_1 @@ -0,0 +1,76 @@ +#!/bin/sh + +cat > partialator_merge_check_1.stream << EOF +CrystFEL stream format 2.1 +Command line: indexamajig -i dummy.lst -o dummy.stream --kraken=prawn +----- Begin chunk ----- +Image filename: dummy.h5 +photon_energy_eV = 2000.0 +--- Begin crystal +Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg +astar = -0.0283891 +0.0149254 -0.0257273 nm^-1 +bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1 +cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1 +Reflections measured after indexing + h k l I phase sigma(I) counts fs/px ss/px + 1 0 0 100.00 - 1.00 1 938.0 629.0 +End of reflections +--- End crystal +--- Begin crystal +Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg +astar = -0.0283891 +0.0149254 -0.0257273 nm^-1 +bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1 +cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1 +Reflections measured after indexing + h k l I phase sigma(I) counts fs/px ss/px + 1 0 0 200.00 - 1.00 1 938.0 629.0 +End of reflections +--- End crystal +----- End chunk ----- +EOF + +cat > partialator_merge_check_1_ans.hkl << EOF + h k l I phase sigma(I) counts fs/px ss/px + 1 0 0 150.00 - 1.00 1 938.0 629.0 +End of reflections +EOF + +cat > partialator_merge_check_1.beam << EOF +beam/fluence = 2.0e15 +beam/radius = 1.5e-6 +beam/photon_energy = 6000.0 +beam/bandwidth = 0.0005 +beam/divergence = 0.001 +profile_radius = 0.005e9 +EOF + +cat > partialator_merge_check_1.geom << EOF +0/min_fs = 0 +0/max_fs = 1023 +0/min_ss = 0 +0/max_ss = 1023 +0/corner_x = -512.00 +0/corner_y = -512.00 +0/fs = x +0/ss = y +0/clen = 70.0e-3 +0/res = 13333.3 ; 75 micron pixel size +0/badrow_direction = y +0/adu_per_eV = 1.0 +EOF + + +src/partialator -i partialator_merge_check_1.stream \ + -o partialator_merge_check_1.hkl \ + -g partialator_merge_check_1.geom \ + -b partialator_merge_check_1.beam \ + --model=unity --iterations=0 --no-scale + +diff partialator_merge_check_1.hkl partialator_merge_check_1_ans.hkl +if [ $? -ne 0 ]; then + exit 1 +fi +rm -f partialator_merge_check_1.stream partialator_merge_check_1.hkl \ + partialator_merge_check_1_ans.hkl partialator_merge_check_1.beam \ + partialator_merge_check_1.geom +exit 0 -- cgit v1.2.3 From f93b8b6546543274e9178f548ef3cc7710390e64 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 30 Apr 2013 17:24:11 +0200 Subject: Fix partialator_merge_check_1 --- src/hrs-scaling.c | 2 +- tests/partialator_merge_check_1 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index 5acb247f..3c0709b4 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -498,7 +498,7 @@ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref, int i; double max_corr; RefList *full = NULL; - const int min_redundancy = 3; + const int min_redundancy = 2; for ( i=0; i partialator_merge_check_1_ans.hkl << EOF h k l I phase sigma(I) counts fs/px ss/px - 1 0 0 150.00 - 1.00 1 938.0 629.0 + 1 0 0 150.00 - 35.36 2 0.0 0.0 End of reflections EOF -- cgit v1.2.3 From 08e0c29ee81bc3602e89b5bf61ac75d200e90239 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 30 Apr 2013 17:48:57 +0200 Subject: Add new tests --- Makefile.am | 13 +++++-- tests/partialator_merge_check_1 | 2 + tests/partialator_merge_check_2 | 79 ++++++++++++++++++++++++++++++++++++++++ tests/partialator_merge_check_3 | 81 +++++++++++++++++++++++++++++++++++++++++ 4 files changed, 171 insertions(+), 4 deletions(-) create mode 100755 tests/partialator_merge_check_2 create mode 100755 tests/partialator_merge_check_3 diff --git a/Makefile.am b/Makefile.am index f465ba36..af810dd8 100644 --- a/Makefile.am +++ b/Makefile.am @@ -11,14 +11,19 @@ noinst_PROGRAMS = tests/list_check tests/integration_check \ tests/centering_check tests/transformation_check \ tests/cell_check -TESTS = tests/list_check tests/first_merge_check tests/second_merge_check \ - tests/third_merge_check tests/fourth_merge_check \ +MERGE_CHECKS = tests/first_merge_check tests/second_merge_check \ + tests/third_merge_check tests/fourth_merge_check + +PARTIAL_CHECKS = tests/partialator_merge_check_1 \ + tests/partialator_merge_check_2 \ + tests/partialator_merge_check_3 + +TESTS = tests/list_check $(MERGE_CHECKS) $(PARTIAL_CHECKS) \ tests/integration_check tests/pr_gradient_check \ tests/symmetry_check tests/centering_check tests/transformation_check \ tests/cell_check -EXTRA_DIST += tests/first_merge_check tests/second_merge_check \ - tests/third_merge_check tests/fourth_merge_check +EXTRA_DIST += $(MERGE_CHECKS) $(PARTIAL_CHECKS) if BUILD_HDFSEE bin_PROGRAMS += src/hdfsee diff --git a/tests/partialator_merge_check_1 b/tests/partialator_merge_check_1 index 98c34d65..9125b75b 100755 --- a/tests/partialator_merge_check_1 +++ b/tests/partialator_merge_check_1 @@ -29,6 +29,8 @@ End of reflections ----- End chunk ----- EOF +# We merge two patterns, without scaling or partiality, the result should just +# be an average. cat > partialator_merge_check_1_ans.hkl << EOF h k l I phase sigma(I) counts fs/px ss/px 1 0 0 150.00 - 35.36 2 0.0 0.0 diff --git a/tests/partialator_merge_check_2 b/tests/partialator_merge_check_2 new file mode 100755 index 00000000..50b0c9bf --- /dev/null +++ b/tests/partialator_merge_check_2 @@ -0,0 +1,79 @@ +#!/bin/sh + +cat > partialator_merge_check_2.stream << EOF +CrystFEL stream format 2.1 +Command line: indexamajig -i dummy.lst -o dummy.stream --kraken=prawn +----- Begin chunk ----- +Image filename: dummy.h5 +photon_energy_eV = 2000.0 +--- Begin crystal +Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg +astar = -0.0283891 +0.0149254 -0.0257273 nm^-1 +bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1 +cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1 +Reflections measured after indexing + h k l I phase sigma(I) counts fs/px ss/px + 1 0 0 100.00 - 1.00 1 938.0 629.0 +End of reflections +--- End crystal +--- Begin crystal +Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg +astar = -0.0283891 +0.0149254 -0.0257273 nm^-1 +bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1 +cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1 +Reflections measured after indexing + h k l I phase sigma(I) counts fs/px ss/px + 1 0 0 200.00 - 1.00 1 938.0 629.0 +End of reflections +--- End crystal +----- End chunk ----- +EOF + +# We merge two patterns, without partiality but with scaling, the result should +# be the mean but with the standard deviation should be zero because the scaling +# factor can absorb the difference. +cat > partialator_merge_check_2_ans.hkl << EOF + h k l I phase sigma(I) counts fs/px ss/px + 1 0 0 150.00 - 0.00 2 0.0 0.0 +End of reflections +EOF + +cat > partialator_merge_check_2.beam << EOF +beam/fluence = 2.0e15 +beam/radius = 1.5e-6 +beam/photon_energy = 6000.0 +beam/bandwidth = 0.0005 +beam/divergence = 0.001 +profile_radius = 0.005e9 +EOF + +cat > partialator_merge_check_2.geom << EOF +0/min_fs = 0 +0/max_fs = 1023 +0/min_ss = 0 +0/max_ss = 1023 +0/corner_x = -512.00 +0/corner_y = -512.00 +0/fs = x +0/ss = y +0/clen = 70.0e-3 +0/res = 13333.3 ; 75 micron pixel size +0/badrow_direction = y +0/adu_per_eV = 1.0 +EOF + + +src/partialator -i partialator_merge_check_2.stream \ + -o partialator_merge_check_2.hkl \ + -g partialator_merge_check_2.geom \ + -b partialator_merge_check_2.beam \ + --model=unity --iterations=0 + +diff partialator_merge_check_2.hkl partialator_merge_check_2_ans.hkl +if [ $? -ne 0 ]; then + exit 1 +fi +rm -f partialator_merge_check_2.stream partialator_merge_check_2.hkl \ + partialator_merge_check_2_ans.hkl partialator_merge_check_2.beam \ + partialator_merge_check_2.geom +exit 0 diff --git a/tests/partialator_merge_check_3 b/tests/partialator_merge_check_3 new file mode 100755 index 00000000..82b3b787 --- /dev/null +++ b/tests/partialator_merge_check_3 @@ -0,0 +1,81 @@ +#!/bin/sh + +cat > partialator_merge_check_3.stream << EOF +CrystFEL stream format 2.1 +Command line: indexamajig -i dummy.lst -o dummy.stream --kraken=prawn +----- Begin chunk ----- +Image filename: dummy.h5 +photon_energy_eV = 2000.0 +--- Begin crystal +Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg +astar = -0.0283891 +0.0149254 -0.0257273 nm^-1 +bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1 +cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1 +Reflections measured after indexing + h k l I phase sigma(I) counts fs/px ss/px + 0 1 0 110.00 - 1.00 1 938.0 629.0 + 2 0 0 190.00 - 1.00 1 938.0 629.0 +End of reflections +--- End crystal +--- Begin crystal +Cell parameters 27.74398 27.84377 16.90346 nm, 88.53688 91.11774 118.75944 deg +astar = -0.0283891 +0.0149254 -0.0257273 nm^-1 +bstar = -0.0068281 +0.0403989 -0.0005196 nm^-1 +cstar = +0.0406926 +0.0052233 -0.0426520 nm^-1 +Reflections measured after indexing + h k l I phase sigma(I) counts fs/px ss/px + 1 0 0 50.00 - 1.00 1 938.0 629.0 + 0 2 0 100.00 - 1.00 1 938.0 629.0 + 0 1 1 1.00 - 1.00 1 938.0 629.0 +End of reflections +--- End crystal +----- End chunk ----- +EOF + +# W +cat > partialator_merge_check_3_ans.hkl << EOF + h k l I phase sigma(I) counts fs/px ss/px + 1 0 0 81.17 - 3.62 2 0.0 0.0 + 2 0 0 144.36 - 2.04 2 0.0 0.0 +End of reflections +EOF + +cat > partialator_merge_check_3.beam << EOF +beam/fluence = 2.0e15 +beam/radius = 1.5e-6 +beam/photon_energy = 6000.0 +beam/bandwidth = 0.0005 +beam/divergence = 0.001 +profile_radius = 0.005e9 +EOF + +cat > partialator_merge_check_3.geom << EOF +0/min_fs = 0 +0/max_fs = 1023 +0/min_ss = 0 +0/max_ss = 1023 +0/corner_x = -512.00 +0/corner_y = -512.00 +0/fs = x +0/ss = y +0/clen = 70.0e-3 +0/res = 13333.3 ; 75 micron pixel size +0/badrow_direction = y +0/adu_per_eV = 1.0 +EOF + + +src/partialator -i partialator_merge_check_3.stream \ + -o partialator_merge_check_3.hkl \ + -g partialator_merge_check_3.geom \ + -b partialator_merge_check_3.beam \ + --model=unity --iterations=0 -y 4 + +diff partialator_merge_check_3.hkl partialator_merge_check_3_ans.hkl +if [ $? -ne 0 ]; then + exit 1 +fi +rm -f partialator_merge_check_3.stream partialator_merge_check_3.hkl \ + partialator_merge_check_3_ans.hkl partialator_merge_check_3.beam \ + partialator_merge_check_3.geom +exit 0 -- cgit v1.2.3 From d2675aac557f43909fde3070ca5ee068ff614fc3 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 30 Apr 2013 18:00:56 +0200 Subject: Remove unused line --- src/hrs-scaling.c | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index 3c0709b4..46590f7e 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -538,8 +538,6 @@ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref, full = lsq_intensities(crystals, n, n_threads, pmodel); } - //show_scale_factors(images, n); - i++; } while ( (max_corr > 0.01) && (i < MAX_CYCLES) ); -- cgit v1.2.3 From 6b18fa18fbdc1f272e65f81d04b22d973acafa1e Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 30 Apr 2013 18:01:07 +0200 Subject: partialator: Add --min-measurements --- src/hrs-scaling.c | 4 ++-- src/hrs-scaling.h | 3 ++- src/partialator.c | 19 +++++++++++++++++-- 3 files changed, 21 insertions(+), 5 deletions(-) diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index 46590f7e..c249ef46 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -493,12 +493,12 @@ static void calculate_esds(Crystal **crystals, int n, RefList *full, /* Scale the stack of images */ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref, - int n_threads, int noscale, PartialityModel pmodel) + int n_threads, int noscale, PartialityModel pmodel, + int min_redundancy) { int i; double max_corr; RefList *full = NULL; - const int min_redundancy = 2; for ( i=0; i Specify partiality model.\n" +" --min-measurements= Require at least measurements before a\n" +" reflection appears in the output. Default: 2\n" "\n" " -j Run analyses in parallel.\n"); } @@ -323,6 +325,8 @@ int main(int argc, char *argv[]) Crystal **crystals; char *pmodel_str = NULL; PartialityModel pmodel = PMODEL_SPHERE; + int min_measurements = 2; + char *rval; /* Long options */ const struct option longopts[] = { @@ -336,6 +340,7 @@ int main(int argc, char *argv[]) {"no-scale", 0, &noscale, 1}, {"reference", 1, NULL, 'r'}, {"model", 1, NULL, 'm'}, + {"min-measurements", 1, NULL, 2}, {0, 0, NULL, 0} }; @@ -398,6 +403,15 @@ int main(int argc, char *argv[]) reference_file = strdup(optarg); break; + case 2 : + errno = 0; + min_measurements = strtod(optarg, &rval); + if ( *rval != '\0' ) { + ERROR("Invalid value for --min-measurements.\n"); + return 1; + } + break; + case 0 : break; @@ -583,7 +597,7 @@ int main(int argc, char *argv[]) STATUS("Performing initial scaling.\n"); if ( noscale ) STATUS("Scale factors fixed at 1.\n"); full = scale_intensities(crystals, n_crystals, reference, - nthreads, noscale, pmodel); + nthreads, noscale, pmodel, min_measurements); sr = sr_titlepage(crystals, n_crystals, "scaling-report.pdf", infile, cmdline); @@ -621,7 +635,8 @@ int main(int argc, char *argv[]) /* Re-estimate all the full intensities */ reflist_free(full); full = scale_intensities(crystals, n_crystals, - reference, nthreads, noscale, pmodel); + reference, nthreads, noscale, pmodel, + min_measurements); sr_iteration(sr, i+1, crystals, n_crystals, full); -- cgit v1.2.3