aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2018-03-29 16:00:50 +0200
committerThomas White <taw@physics.org>2018-03-29 16:00:50 +0200
commit8d9f9cc570a69780d3e3b2c271f03d055066a9be (patch)
tree44d7a71e0cbcab2c1c327d761da150bca9c1807a
parent77cf2edd09bb01ae331935f467064c751f6e338e (diff)
parent402b3870e3e1bceb974ad5a402c5e5e898f4c87e (diff)
Merge branch 'tom/ginn-partials'
Diffstat
-rw-r--r--Makefile.am18
-rw-r--r--doc/man/partialator.146
-rw-r--r--libcrystfel/src/geometry.c465
-rw-r--r--libcrystfel/src/geometry.h21
-rw-r--r--libcrystfel/src/image.h2
-rw-r--r--libcrystfel/src/integration.c2
-rw-r--r--libcrystfel/src/predict-refine.c8
-rw-r--r--libcrystfel/src/reflist.c92
-rw-r--r--libcrystfel/src/reflist.h13
-rwxr-xr-xscripts/plot-contourmap153
-rwxr-xr-xscripts/plot-pr196
-rw-r--r--src/merge.c13
-rw-r--r--src/merge.h3
-rw-r--r--src/partial_sim.c5
-rw-r--r--src/partialator.c448
-rw-r--r--src/post-refinement.c1178
-rw-r--r--src/post-refinement.h21
-rw-r--r--src/process_image.c4
-rw-r--r--src/rejection.c2
-rw-r--r--src/scaling.c215
-rw-r--r--src/scaling.h13
-rw-r--r--tests/.gitignore1
-rw-r--r--tests/linear_scale_check.c77
-rw-r--r--tests/pr_p_gradient_check.c529
-rw-r--r--tests/prediction_gradient_check.c22
25 files changed, 2136 insertions, 1411 deletions
diff --git a/Makefile.am b/Makefile.am
index c3072602..94ec2d6f 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -8,25 +8,25 @@ bin_PROGRAMS = src/pattern_sim src/process_hkl src/get_hkl src/indexamajig \
src/ambigator src/geoptimiser src/whirligig src/list_events
noinst_PROGRAMS = tests/list_check tests/integration_check \
- tests/pr_p_gradient_check tests/symmetry_check \
+ tests/symmetry_check \
tests/centering_check tests/transformation_check \
tests/cell_check tests/ring_check \
tests/prof2d_check tests/ambi_check \
- tests/prediction_gradient_check
+ tests/prediction_gradient_check \
+ tests/linear_scale_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/pr_p_gradient_check
+ tests/partialator_merge_check_3
TESTS = tests/list_check $(MERGE_CHECKS) $(PARTIAL_CHECKS) \
tests/integration_check \
tests/symmetry_check tests/centering_check tests/transformation_check \
tests/cell_check tests/ring_check tests/prof2d_check tests/ambi_check \
- tests/prediction_gradient_check
+ tests/prediction_gradient_check tests/linear_scale_check
EXTRA_DIST += $(MERGE_CHECKS) $(PARTIAL_CHECKS)
EXTRA_DIST += relnotes-0.6.3 announcement-0.6.3
@@ -113,9 +113,6 @@ tests_symmetry_check_SOURCES = tests/symmetry_check.c
tests_ambi_check_SOURCES = tests/ambi_check.c
-tests_pr_p_gradient_check_SOURCES = tests/pr_p_gradient_check.c \
- src/post-refinement.c
-
tests_prediction_gradient_check_SOURCES = tests/prediction_gradient_check.c
tests_centering_check_SOURCES = tests/centering_check.c
@@ -126,6 +123,9 @@ tests_ring_check_SOURCES = tests/ring_check.c
tests_cell_check_SOURCES = tests/cell_check.c
+tests_linear_scale_check_SOURCES = tests/linear_scale_check.c src/scaling.c \
+ src/merge.c
+
INCLUDES = -I$(top_srcdir)/libcrystfel/src -I$(top_srcdir)/data
EXTRA_DIST += src/dw-hdfsee.h src/hdfsee.h src/render_hkl.h \
@@ -185,7 +185,7 @@ script_DATA = scripts/alternate-stream scripts/cell-please \
scripts/gaincal-to-saturation-map scripts/move-entire-detector \
scripts/split-by-mask scripts/turbo-index-slurm \
scripts/sum-peaks scripts/peakogram-stream scripts/eiger-badmap \
- scripts/sum-hdf5-files
+ scripts/sum-hdf5-files scripts/plot-pr scripts/plot-contourmap
EXTRA_DIST += $(script_DATA)
diff --git a/doc/man/partialator.1 b/doc/man/partialator.1
index decaa07c..0a5b78ba 100644
--- a/doc/man/partialator.1
+++ b/doc/man/partialator.1
@@ -1,7 +1,7 @@
.\"
.\" partialator man page
.\"
-.\" Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
+.\" Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
.\" a research centre of the Helmholtz Association.
.\"
.\" Part of CrystFEL - crystallography with a FEL
@@ -132,29 +132,45 @@ Reject crystals if the absolute values of their relative Debye-Waller ("B") fact
.PD
Write out the per-crystal parameters and reflection lists after every cycle of refinement, instead of only at the end. The intermediate reflection lists and parameter filenames will be prefixed with \fBiter\fIN\fB_\fR (note the trailing underscore), where \fIN\fR is the iteration number. If you use \fB--custom-split\fR, intermediate results will also be output for each custom dataset.
+.PD 0
+.IP \fB--no-logs\fR
+.PD
+Do not write the extensive log files needed for plotting contour maps and spectrum graphs. This makes the process a lot faster, but you probably do want these logs to check that post-refinement is working reasonably.
+
+.PD 0
+.IP "\fB-w\fR \fIpg\fR"
+.PD
+Set the apparent symmetry of the crystals. The ambiguity operator will be determined by comparing this to the actual symmetry.
+.IP
+If you prefer (or the scenario demands it), you can specify the ambiguity operator directly using \fB--operator\fR.
+
+.PD 0
+.IP \fB--operator=\fR\fIop\fR
+.PD
+Specify the indexing ambiguity operator. Example: \fB--operator=k,h,-l\fR.
+.IP
+If you prefer, you can specify the ambiguity operator by specifying the apparent symmetry using \fB-w\fR.
+
+.PD 0
+.IP \fB--force-bandwidth=\fIbw\fR
+.IP \fB--force-radius=\fIR\fR
+.PD
+Set the X-ray bandwidth or initial profile radius for all crystals before proceeding. Bandwidth is given as a fraction, i.e. \fB--force-bandwidth=0.0013\fR means 0.13 percent (approximate FWHM). Radius is given in nm^-1.
+
.SH PARTIALITY MODELS
The available partiality models are:
-.IP \fBscsphere\fR
+.IP \fBxsphere\fR
.PD
The volume of intersection between a sphere centered on each reciprocal lattice
point and the part of reciprocal space excited by the Ewald sphere taking into
account the finite bandwidth and convergence angle. A "source coverage factor"
is included to take into account the spectral brightness of the effective
-source for the reflection.
-
-This model is similar to that described in Acta Cryst. D69 (2013) p1231-1240,
-and in Phil. Trans. Roy. Soc. B 369 (2014) 20130330, except that the "Lorentz
-factor" described there is no longer treated as a separate factor.
+source for the reflection. The X-ray spectrum is modelled as a super-Gaussian
+with exponent 1.5, and the overlap integral is evaluated numerically.
-
-.IP \fBscgaussian\fR
-.PD
-As \fBscsphere\fR, except that the shape of the scattering density centered on
-each reciprocal lattice point is taken to be a 3D Gaussian distribution instead
-of a sphere. The standard deviation of the distribution will be the profile
-radius (determined by indexamajig) divided by 2.6.
+This model is the same as that described in Acta Cryst. D71 (2015) p1400.
.IP \fBunity\fR
.PD
@@ -216,7 +232,7 @@ This page was written by Thomas White.
Report bugs to <taw@physics.org>, or visit <http://www.desy.de/~twhite/crystfel>.
.SH COPYRIGHT AND DISCLAIMER
-Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, a research centre of the Helmholtz Association.
+Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, a research centre of the Helmholtz Association.
.P
partialator, and this manual, are part of CrystFEL.
.P
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 576655de..9ae6d96f 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -238,76 +238,38 @@ static double random_partiality(signed int h, signed int k, signed int l,
}
-static double partiality(PartialityModel pmodel,
- signed int h, signed int k, signed int l,
- int serial,
- double rlow, double rhigh, double pr)
-{
- double D = rlow - rhigh;
-
- /* Convert to partiality */
- switch ( pmodel ) {
-
- default:
- case PMODEL_UNITY:
- return 1.0;
-
- case PMODEL_SCSPHERE:
- return 4.0*sphere_fraction(rlow, rhigh, pr)*pr / (3.0*D);
-
- case PMODEL_SCGAUSSIAN:
- return 4.0*gaussian_fraction(rlow, rhigh, pr)*pr / (3.0*D);
-
- case PMODEL_RANDOM:
- return random_partiality(h, k, l, serial);
-
- }
-}
-
-
static Reflection *check_reflection(struct image *image, Crystal *cryst,
signed int h, signed int k, signed int l,
double xl, double yl, double zl,
Reflection *updateme)
{
- double tl;
- double rlow, rhigh; /* "Excitation error" */
- double klow, khigh; /* Wavenumber */
Reflection *refl;
- double cet, cez; /* Centre of Ewald sphere */
- double pr;
- double del;
+ double R, top;
+ double kmin, kmax, k0, knom, k1, khalf;
+ double dcs, exerr;
/* Don't predict 000 */
if ( (updateme == NULL) && (abs(h)+abs(k)+abs(l) == 0) ) return NULL;
- pr = crystal_get_profile_radius(cryst);
- del = image->div + crystal_get_mosaicity(cryst);
-
- /* "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);
-
- /* If the point is looking "backscattery", reject it straight away */
- if ( (updateme == NULL) && (zl < -khigh/2.0) ) return NULL;
-
- tl = sqrt(xl*xl + yl*yl);
-
- cet = -sin(del/2.0) * khigh;
- cez = -cos(del/2.0) * khigh;
- rhigh = khigh - distance(cet, cez, tl, zl); /* Loss of precision */
-
- cet = sin(del/2.0) * klow;
- cez = -cos(del/2.0) * klow;
- rlow = klow - distance(cet, cez, tl, zl); /* Loss of precision */
-
- /* Condition for reflection to be excited at all */
- if ( (updateme == NULL)
- && (signbit(rlow) == signbit(rhigh))
- && (fabs(rlow) > pr)
- && (fabs(rhigh) > pr) ) return NULL;
+ /* Calculate the limiting wavelengths, lambda0 and lambda1
+ * = 1/k0 and 1/k1 respectively */
+ R = fabs(crystal_get_profile_radius(cryst));
+ top = R*R - xl*xl - yl*yl - zl*zl;
+ k0 = top/(2.0*(zl+R));
+ khalf = (- xl*xl - yl*yl - zl*zl) / (2.0*zl);
+ k1 = top/(2.0*(zl-R));
+
+ /* The reflection is excited if any of the reflection is within 2sigma
+ * of the nominal * wavelength of the X-ray beam
+ * (NB image->bw is full width) */
+ kmin = 1.0/(image->lambda + image->lambda*image->bw);
+ knom = 1.0/image->lambda;
+ kmax = 1.0/(image->lambda - image->lambda*image->bw);
+ if ( (updateme == NULL) && ((k1>kmax) || (k0<kmin)) ) return NULL;
+
+ /* Calculate excitation error */
+ dcs = distance3d(0.0, 0.0, -knom, xl, yl, zl);
+ exerr = 1.0/image->lambda - dcs; /* Loss of precision */
if ( updateme == NULL ) {
refl = reflection_new(h, k, l);
@@ -321,7 +283,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
if ( (image->det != NULL) && (updateme != NULL) ) {
double fs, ss;
- locate_peak_on_panel(xl, yl, zl, 1.0/image->lambda,
+ locate_peak_on_panel(xl, yl, zl, knom,
get_panel(updateme), &fs, &ss);
set_detector_pos(refl, fs, ss);
@@ -333,7 +295,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
double fs, ss; /* Position on detector */
signed int p; /* Panel number */
- p = locate_peak(xl, yl, zl, 1.0/image->lambda,
+ p = locate_peak(xl, yl, zl, knom,
image->det, &fs, &ss);
if ( p == -1 ) {
reflection_free(refl);
@@ -344,20 +306,9 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
}
- if ( unlikely(rlow < rhigh) ) {
- ERROR("Reflection with rlow < rhigh!\n");
- ERROR("%3i %3i %3i rlow = %e, rhigh = %e\n",
- h, k, l, rlow, rhigh);
- ERROR("div + m = %e, R = %e, bw = %e\n", del, pr, image->bw);
- /* If we are updating, this is (kind of) OK */
- if ( updateme == NULL ) {
- reflection_free(refl);
- return NULL;
- }
- }
-
- set_partial(refl, rlow, rhigh, 1.0); /* Actual partiality set by
- * calculate_partialities() */
+ set_kpred(refl, knom);
+ set_khalf(refl, khalf);
+ set_exerr(refl, exerr);
set_lorentz(refl, 1.0);
set_symmetric_indices(refl, h, k, l);
set_redundancy(refl, 1);
@@ -368,18 +319,12 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
double r_gradient(UnitCell *cell, int k, Reflection *refl, struct image *image)
{
- double azi;
double asx, asy, asz;
double bsx, bsy, bsz;
double csx, csy, csz;
double xl, yl, zl;
signed int hs, ks, ls;
- double rlow, rhigh, p;
- double philow, phihigh, phi;
- double khigh, klow;
- double tl, cet, cez;
-
- get_partial(refl, &rlow, &rhigh, &p);
+ double tl, phi, azi;
get_symmetric_indices(refl, &hs, &ks, &ls);
@@ -390,26 +335,9 @@ double r_gradient(UnitCell *cell, int k, Reflection *refl, struct image *image)
yl = hs*asy + ks*bsy + ls*csy;
zl = hs*asz + ks*bsz + ls*csz;
- /* "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);
-
tl = sqrt(xl*xl + yl*yl);
-
- cet = -sin(image->div/2.0) * klow;
- cez = -cos(image->div/2.0) * klow;
- philow = angle_between_2d(tl-cet, zl-cez, 0.0, 1.0);
-
- cet = -sin(image->div/2.0) * khigh;
- cez = -cos(image->div/2.0) * khigh;
- phihigh = angle_between_2d(tl-cet, zl-cez, 0.0, 1.0);
-
- /* Approximation: philow and phihigh are very similar */
- phi = (philow + phihigh) / 2.0;
-
- azi = atan2(yl, xl);
+ phi = angle_between_2d(tl, zl+1.0/image->lambda, 0.0, 1.0); /* 2theta */
+ azi = atan2(yl, xl); /* azimuth */
switch ( k ) {
@@ -531,11 +459,17 @@ RefList *predict_to_res(Crystal *cryst, double max_res)
}
-static void set_unity_partialities(RefList *list)
+static void set_unity_partialities(Crystal *cryst)
{
+ RefList *list;
Reflection *refl;
RefListIterator *iter;
+ list = crystal_get_reflections(cryst);
+ if ( list == NULL ) {
+ ERROR("No reflections for partiality calculation!\n");
+ return;
+ }
for ( refl = first_refl(list, &iter);
refl != NULL;
refl = next_refl(refl, iter) )
@@ -546,25 +480,11 @@ static void set_unity_partialities(RefList *list)
}
-/**
- * calculate_partialities:
- * @cryst: A %Crystal
- * @pmodel: A %PartialityModel
- *
- * Calculates the partialities for the reflections in @cryst, given the current
- * crystal and image parameters. The crystal's image and reflection lists
- * must be set. The specified %PartialityModel will be used.
- *
- * You must not have changed the crystal or image parameters since you last
- * called predict_to_res() or update_predictions(), because this function
- * relies on the limiting wavelength values calculated by those functions.
- */
-void calculate_partialities(Crystal *cryst, PartialityModel pmodel)
+static void set_random_partialities(Crystal *cryst)
{
RefList *list;
Reflection *refl;
RefListIterator *iter;
- double pr;
struct image *image;
list = crystal_get_reflections(cryst);
@@ -573,8 +493,166 @@ void calculate_partialities(Crystal *cryst, PartialityModel pmodel)
return;
}
- if ( pmodel == PMODEL_UNITY ) {
- set_unity_partialities(list);
+ image = crystal_get_image(cryst);
+ if ( image == NULL ) {
+ ERROR("No image structure for partiality calculation!\n");
+ return;
+ }
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ signed int h, k, l;
+ get_symmetric_indices(refl, &h, &k, &l);
+ set_partiality(refl, random_partiality(h, k, l, image->serial));
+ set_lorentz(refl, 1.0);
+ }
+}
+
+
+static double do_integral(double q2, double zl, double R,
+ double lambda, double sig, char *verbose)
+{
+ int i;
+ double kmin, kmax, kstart, kfinis;
+ double inc;
+ double total = 0.0;
+ double k0, k1;
+ const int SAMPLES = 50; /* Number of samples for integration */
+ const double N = 1.5; /* Pointiness of spectrum */
+ FILE *fh = NULL;
+
+ assert(R*R < q2);
+ assert(R > 0.0);
+
+ /* Range over which P is different from zero */
+ k0 = (R*R - q2)/(2.0*(zl+R));
+ k1 = (R*R - q2)/(2.0*(zl-R));
+
+ /* Check for reflections partially "behind the beam" */
+ if ( k0 < 0.0 ) k0 = +INFINITY;
+ if ( k1 < 0.0 ) k1 = +INFINITY;
+
+ /* Range over which E is significantly different from zero */
+ kmin = 1.0 / (lambda + 3.0*sig);
+ kmax = 1.0 / (lambda - 3.0*sig);
+
+ /* Calculate range over which E*P is different from zero */
+ if ( k0 < k1 ) {
+ STATUS("%e %e\n", k0, k1);
+ STATUS("%e %e %e\n", q2, zl, R);
+ STATUS("%e %e\n", kmin, kmax);
+ }
+ assert((isinf(k0) && isinf(k1)) || (k0 > k1));
+ assert(kmax > kmin);
+ if ( kmin < k1 ) {
+ if ( kmax < k1 ) {
+ /* Case 1 */
+ kstart = k1; kfinis = k0;
+ return 0.0;
+ } else if ( kmax < k0 ) {
+ /* Case 2 (kmax > k1)*/
+ kstart = k1; kfinis = kmax;
+ } else {
+ /* Case 3 (kmax > k1 and kmax > k0) */
+ kstart = k1; kfinis = k0;
+ }
+ } else if ( kmin < k0 ) { /* kmin > k1 */
+ if ( kmax < k0 ) {
+ /* Case 4 */
+ kstart = kmin; kfinis = kmax;
+ } else {
+ /* Case 5 (kmax > k0) */
+ kstart = kmin; kfinis = k0;
+ }
+ } else {
+ /* Case 6 (kmin > k1 and (kmax>)kmin > k0) */
+ kstart = k1; kfinis = k0;
+ return 0.0;
+ }
+
+ if ( kstart < 0.0 ) kstart = kmin;
+ if ( kfinis < 0.0 ) kfinis = kmax;
+
+ inc = (kfinis - kstart) / SAMPLES;
+
+ if ( verbose ) {
+ char fn[64];
+ snprintf(fn, 63, "partial%s.graph", verbose);
+ fh = fopen(fn, "w");
+ fprintf(fh, " n p wavelength E P\n");
+ STATUS("Nominal k = %e m^-1\n", 1.0/lambda);
+ STATUS(" (wavelength %e m)\n", lambda);
+ STATUS("Bandwidth %e m\n", sig);
+ STATUS("k1/2 = %e m^-1\n", -q2/(2.0*zl));
+ STATUS(" (wavelength %e m)\n", 1.0/(-q2/(2.0*zl)));
+ STATUS("Reflection k goes from %e to %e m^-1\n", k1, k0);
+ STATUS(" (wavelengths from %e to %e m\n", 1.0/k1, 1.0/k0);
+ STATUS("Beam goes from %e to %e m^-1\n", kmin, kmax);
+ STATUS(" (wavelengths from %e to %e m\n", 1.0/kmin, 1.0/kmax);
+ STATUS("Integration goes from %e to %e m^-1\n", kstart, kfinis);
+ STATUS(" (wavelengths from %e to %e m\n", 1.0/kstart, 1.0/kfinis);
+ }
+
+ for ( i=0; i<SAMPLES; i++ ) {
+
+ double p, kp, lrel;
+ double E, P;
+
+ kp = kstart + i*inc;
+ double pref = sqrt(q2 + kp*kp + 2.0*zl*kp)/(2.0*R);
+ p = pref + 0.5 - kp/(2.0*R);
+
+ /* Spectral energy term */
+ lrel = fabs(1.0/kp - lambda);
+ E = exp(-0.5 * pow(lrel / sig, N));
+ E /= sqrt(2.0 * M_PI * sig);
+
+ /* RLP profile term */
+ P = 4.0*p * (1.0 - p);
+
+ total += E*P*inc;
+
+ if ( fh != NULL ) {
+ fprintf(fh, "%3i %f %e %e %e\n", i, p, 1.0/kp, E, P);
+ }
+ }
+
+ if ( fh != NULL ) fclose(fh);
+
+ if ( isnan(total) ) {
+ ERROR("NaN total!\n");
+ STATUS("Nominal k = %e m^-1\n", 1.0/lambda);
+ STATUS(" (wavelength %e m)\n", lambda);
+ STATUS("Bandwidth %e m\n", sig);
+ STATUS("k1/2 = %e m^-1\n", -q2/(2.0*zl));
+ STATUS(" (wavelength %e m)\n", 1.0/(-q2/(2.0*zl)));
+ STATUS("Reflection k goes from %e to %e m^-1\n", k1, k0);
+ STATUS(" (wavelengths from %e to %e m\n", 1.0/k1, 1.0/k0);
+ STATUS("Beam goes from %e to %e m^-1\n", kmin, kmax);
+ STATUS(" (wavelengths from %e to %e m\n", 1.0/kmin, 1.0/kmax);
+ STATUS("Integration goes from %e to %e m^-1\n", kstart, kfinis);
+ STATUS(" (wavelengths from %e to %e m\n", 1.0/kstart, 1.0/kfinis);
+ }
+
+ return total;
+}
+
+
+static void ginn_spectrum_partialities(Crystal *cryst)
+{
+ RefList *list;
+ Reflection *refl;
+ RefListIterator *iter;
+ double r0, m, lambda, sig;
+ struct image *image;
+ UnitCell *cell;
+ double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
+
+ list = crystal_get_reflections(cryst);
+ if ( list == NULL ) {
+ ERROR("No reflections for partiality calculation!\n");
return;
}
@@ -584,31 +662,88 @@ void calculate_partialities(Crystal *cryst, PartialityModel pmodel)
return;
}
- pr = crystal_get_profile_radius(cryst);
+ cell = crystal_get_cell(cryst);
+ if ( cell == NULL ) {
+ ERROR("No unit cell for partiality calculation!\n");
+ return;
+ }
+ cell_get_reciprocal(cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+
+ r0 = fabs(crystal_get_profile_radius(cryst));
+ m = crystal_get_mosaicity(cryst);
+ lambda = image->lambda;
+ sig = image->bw * lambda / 2.0;
for ( refl = first_refl(crystal_get_reflections(cryst), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
- double rlow, rhigh, part;
+ double R;
signed int h, k, l;
+ double xl, yl, zl;
+ double q2;
+ double total, norm;
get_symmetric_indices(refl, &h, &k, &l);
- get_partial(refl, &rlow, &rhigh, &part);
+ xl = h*asx + k*bsx + l*csx;
+ yl = h*asy + k*bsy + l*csy;
+ zl = h*asz + k*bsz + l*csz;
+
+ /* Radius of rlp profile */
+ q2 = xl*xl + yl*yl + zl*zl;
+
+ R = r0 + m * sqrt(q2);
+
+ //char tmp[256];
+ //snprintf(tmp, 255, "-%i,%i,%i", h, k, l);
+ char *tmp = NULL;
- part = partiality(pmodel, h, k, l, image->serial,
- rlow, rhigh, pr);
+ total = do_integral(q2, zl, R, lambda, sig, tmp);
+ norm = do_integral(q2, -0.5*q2*lambda, R, lambda, sig, NULL);
- if ( isnan(part) && ((h!=0) || (k!=0) || (l!=0)) ) {
- ERROR("Assigning NAN partiality!\n");
- ERROR("%3i %3i %3i rlow = %e, rhigh = %e\n",
- h, k, l, rlow, rhigh);
- ERROR("R = %e, bw = %e\n", pr, image->bw);
- ERROR("D = %e\n", rlow - rhigh);
- abort();
- }
+ set_partiality(refl, total/norm);
+ set_lorentz(refl, 1.0);
- set_partiality(refl, part);
+ assert(total <= 2.0*norm);
+
+ }
+}
+
+
+/**
+ * calculate_partialities:
+ * @cryst: A %Crystal
+ * @pmodel: A %PartialityModel
+ *
+ * Calculates the partialities for the reflections in @cryst, given the current
+ * crystal and image parameters. The crystal's image and reflection lists
+ * must be set. The specified %PartialityModel will be used.
+ *
+ * You must not have changed the crystal or image parameters since you last
+ * called predict_to_res() or update_predictions(), because this function
+ * relies on the limiting wavelength values calculated by those functions.
+ */
+void calculate_partialities(Crystal *cryst, PartialityModel pmodel)
+{
+ switch ( pmodel ) {
+
+ case PMODEL_UNITY :
+ set_unity_partialities(cryst);
+ break;
+
+ case PMODEL_XSPHERE :
+ ginn_spectrum_partialities(cryst);
+ break;
+
+ case PMODEL_RANDOM :
+ set_random_partialities(cryst);
+ break;
+
+ default :
+ ERROR("Unknown partiality model %i\n", pmodel);
+ break;
}
}
@@ -694,29 +829,30 @@ void polarisation_correction(RefList *list, UnitCell *cell, struct image *image)
/* Returns dx_h/dP, where P = any parameter */
-double x_gradient(int param, Reflection *refl, UnitCell *cell,
- struct panel *p, double lambda)
+double x_gradient(int param, Reflection *refl, UnitCell *cell, struct panel *p)
{
signed int h, k, l;
- double x, z, wn;
- double ax, ay, az, bx, by, bz, cx, cy, cz;
+ double xl, zl, kpred;
+ double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
get_indices(refl, &h, &k, &l);
- wn = 1.0 / lambda;
- cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
- x = h*ax + k*bx + l*cx;
- z = h*az + k*bz + l*cz;
+ kpred = get_kpred(refl);
+ cell_get_reciprocal(cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+ xl = h*asx + k*bsx + l*csx;
+ zl = h*asz + k*bsz + l*csz;
switch ( param ) {
case GPARAM_ASX :
- return h * p->clen / (wn+z);
+ return h * p->clen / (kpred + zl);
case GPARAM_BSX :
- return k * p->clen / (wn+z);
+ return k * p->clen / (kpred + zl);
case GPARAM_CSX :
- return l * p->clen / (wn+z);
+ return l * p->clen / (kpred + zl);
case GPARAM_ASY :
return 0.0;
@@ -728,13 +864,13 @@ double x_gradient(int param, Reflection *refl, UnitCell *cell,
return 0.0;
case GPARAM_ASZ :
- return -h * x * p->clen / (wn*wn + 2*wn*z + z*z);
+ return -h * xl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl);
case GPARAM_BSZ :
- return -k * x * p->clen / (wn*wn + 2*wn*z + z*z);
+ return -k * xl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl);
case GPARAM_CSZ :
- return -l * x * p->clen / (wn*wn + 2*wn*z + z*z);
+ return -l * xl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl);
case GPARAM_DETX :
return -1;
@@ -743,7 +879,7 @@ double x_gradient(int param, Reflection *refl, UnitCell *cell,
return 0;
case GPARAM_CLEN :
- return x / (wn+z);
+ return xl / (kpred+zl);
}
@@ -753,18 +889,19 @@ double x_gradient(int param, Reflection *refl, UnitCell *cell,
/* Returns dy_h/dP, where P = any parameter */
-double y_gradient(int param, Reflection *refl, UnitCell *cell,
- struct panel *p, double lambda)
+double y_gradient(int param, Reflection *refl, UnitCell *cell, struct panel *p)
{
signed int h, k, l;
- double y, z, wn;
- double ax, ay, az, bx, by, bz, cx, cy, cz;
+ double yl, zl, kpred;
+ double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
get_indices(refl, &h, &k, &l);
- wn = 1.0 / lambda;
- cell_get_reciprocal(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz);
- y = h*ay + k*by + l*cy;
- z = h*az + k*bz + l*cz;
+ kpred = get_kpred(refl);
+ cell_get_reciprocal(cell, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+ yl = h*asy + k*bsy + l*csy;
+ zl = h*asz + k*bsz + l*csz;
switch ( param ) {
@@ -778,22 +915,22 @@ double y_gradient(int param, Reflection *refl, UnitCell *cell,
return 0.0;
case GPARAM_ASY :
- return h * p->clen / (wn+z);
+ return h * p->clen / (kpred + zl);
case GPARAM_BSY :
- return k * p->clen / (wn+z);
+ return k * p->clen / (kpred + zl);
case GPARAM_CSY :
- return l * p->clen / (wn+z);
+ return l * p->clen / (kpred + zl);
case GPARAM_ASZ :
- return -h * y * p->clen / (wn*wn + 2*wn*z + z*z);
+ return -h * yl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl);
case GPARAM_BSZ :
- return -k * y * p->clen / (wn*wn + 2*wn*z + z*z);
+ return -k * yl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl);
case GPARAM_CSZ :
- return -l * y * p->clen / (wn*wn + 2*wn*z + z*z);
+ return -l * yl * p->clen / (kpred*kpred + 2.0*kpred*zl + zl*zl);
case GPARAM_DETX :
return 0;
@@ -802,7 +939,7 @@ double y_gradient(int param, Reflection *refl, UnitCell *cell,
return -1;
case GPARAM_CLEN :
- return y / (wn+z);
+ return yl / (kpred+zl);
}
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index 4fbd0bca..529f112e 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -3,7 +3,7 @@
*
* Geometry of diffraction
*
- * Copyright © 2013-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2013-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
* Copyright © 2012 Richard Kirian
*
@@ -47,8 +47,7 @@ extern "C" {
/**
* PartialityModel:
* @PMODEL_UNITY : Set all all partialities and Lorentz factors to 1.
- * @PMODEL_SCSPHERE : Sphere model with source coverage factor included
- * @PMODEL_SCGAUSSIAN : Gaussian model with source coverage factor included
+ * @PMODEL_XSPHERE : Flat sphere model with super-Gaussian spectrum
* @PMODEL_RANDOM : Randomly assigned partialities
*
* A %PartialityModel describes a geometrical model which can be used to
@@ -57,8 +56,7 @@ extern "C" {
typedef enum {
PMODEL_UNITY,
- PMODEL_SCSPHERE,
- PMODEL_SCGAUSSIAN,
+ PMODEL_XSPHERE,
PMODEL_RANDOM,
} PartialityModel;
@@ -80,8 +78,13 @@ enum gparam {
GPARAM_DETX,
GPARAM_DETY,
GPARAM_CLEN,
- GPARAM_OSF,
- GPARAM_BFAC
+ GPARAM_OSF, /* Linear scale factor */
+ GPARAM_BFAC, /* D-W scale factor */
+ GPARAM_ANG1, /* Out of plane rotation angles of crystal */
+ GPARAM_ANG2,
+ GPARAM_WAVELENGTH,
+
+ GPARAM_EOL /* End of list */
};
@@ -99,9 +102,9 @@ extern double sphere_fraction(double rlow, double rhigh, double pr);
extern double gaussian_fraction(double rlow, double rhigh, double pr);
extern double x_gradient(int param, Reflection *refl, UnitCell *cell,
- struct panel *p, double lambda);
+ struct panel *p);
extern double y_gradient(int param, Reflection *refl, UnitCell *cell,
- struct panel *p, double lambda);
+ struct panel *p);
#ifdef __cplusplus
}
diff --git a/libcrystfel/src/image.h b/libcrystfel/src/image.h
index 99c3d57e..9769c66a 100644
--- a/libcrystfel/src/image.h
+++ b/libcrystfel/src/image.h
@@ -216,7 +216,7 @@ struct image {
/* Per-shot radiation values */
double lambda; /* Wavelength in m */
double div; /* Divergence in radians */
- double bw; /* Bandwidth as a fraction */
+ double bw; /* FWHM bandwidth as a fraction */
/* Detected peaks */
long long num_peaks;
diff --git a/libcrystfel/src/integration.c b/libcrystfel/src/integration.c
index c18ae110..fcbdf658 100644
--- a/libcrystfel/src/integration.c
+++ b/libcrystfel/src/integration.c
@@ -1760,7 +1760,7 @@ void integrate_all_2(struct image *image, IntegrationMethod meth,
IntDiag int_diag,
signed int idh, signed int idk, signed int idl)
{
- integrate_all_3(image, meth, PMODEL_SCSPHERE, 0.0,
+ integrate_all_3(image, meth, PMODEL_XSPHERE, 0.0,
ir_inn, ir_mid, ir_out,
int_diag, idh, idk, idl);
}
diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c
index 3246dbec..14a60bc7 100644
--- a/libcrystfel/src/predict-refine.c
+++ b/libcrystfel/src/predict-refine.c
@@ -91,9 +91,7 @@ static void twod_mapping(double fs, double ss, double *px, double *py,
static double r_dev(struct reflpeak *rp)
{
/* Excitation error term */
- double rlow, rhigh, p;
- get_partial(rp->refl, &rlow, &rhigh, &p);
- return (rlow+rhigh)/2.0;
+ return get_exerr(rp->refl);
}
@@ -425,7 +423,7 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell,
/* Positional x terms */
for ( k=0; k<num_params; k++ ) {
gradients[k] = x_gradient(rv[k], rps[i].refl, cell,
- rps[i].panel, image->lambda);
+ rps[i].panel);
}
for ( k=0; k<num_params; k++ ) {
@@ -457,7 +455,7 @@ static int iterate(struct reflpeak *rps, int n, UnitCell *cell,
/* Positional y terms */
for ( k=0; k<num_params; k++ ) {
gradients[k] = y_gradient(rv[k], rps[i].refl, cell,
- rps[i].panel, image->lambda);
+ rps[i].panel);
}
for ( k=0; k<num_params; k++ ) {
diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c
index 707c802e..a35aa575 100644
--- a/libcrystfel/src/reflist.c
+++ b/libcrystfel/src/reflist.c
@@ -67,10 +67,11 @@ struct _refldata {
signed int ls;
/* Partiality and related geometrical stuff */
- double rlow; /* Low excitation error */
- double rhigh; /* High excitation error */
- double p; /* Partiality */
- double L; /* Lorentz factor */
+ double khalf; /* Wavenumber of middle of reflection */
+ double kpred; /* Wavenumber for prediction */
+ double exerr; /* Excitation error */
+ double p; /* Partiality */
+ double L; /* Lorentz factor */
/* Location in image */
double fs;
@@ -422,22 +423,43 @@ double get_intensity(const Reflection *refl)
/**
- * get_partial:
+ * get_khalf
* @refl: A %Reflection
- * @rlow: Location at which to store the "low" excitation error
- * @rhigh: Location at which to store the "high" excitation error
- * @p: Location at which to store the partiality
*
- * This function is used during post refinement (in conjunction with
- * set_partial()) to get access to the details of the partiality calculation.
+ * Returns: the wavenumber at the centre of the reflection
*
**/
-void get_partial(const Reflection *refl, double *rlow, double *rhigh,
- double *p)
+double get_khalf(const Reflection *refl)
{
- *rlow = refl->data.rlow;
- *rhigh = refl->data.rhigh;
- *p = get_partiality(refl);
+ return refl->data.khalf;
+}
+
+
+
+
+/**
+ * get_kpred:
+ * @refl: A %Reflection
+ *
+ * Returns: the wavenumber which should be used for prediction of this reflection
+ *
+ **/
+double get_kpred(const Reflection *refl)
+{
+ return refl->data.kpred;
+}
+
+
+/**
+ * get_exerr:
+ * @refl: A %Reflection
+ *
+ * Returns: the excitation error (in m^-1) for this reflection
+ *
+ **/
+double get_exerr(const Reflection *refl)
+{
+ return refl->data.exerr;
}
@@ -612,23 +634,43 @@ void set_panel(Reflection *refl, struct panel *p)
refl->data.panel = p;
}
+/**
+ * set_khalf:
+ * @refl: A %Reflection
+ * @khalf: The wavenumber at which the reflection should be predicted
+ *
+ * Sets the wavenumber at the centre of the reflection.
+ **/
+void set_khalf(Reflection *refl, double khalf)
+{
+ refl->data.khalf = khalf;
+}
+
+
/**
- * set_partial:
+ * set_kpred:
* @refl: A %Reflection
- * @rlow: The "low" excitation error
- * @rhigh: The "high" excitation error
- * @p: The partiality
+ * @kpred: The wavenumber at which the reflection should be predicted
*
- * This function is used during post refinement (in conjunction with
- * get_partial()) to get access to the details of the partiality calculation.
+ * Sets the wavenumber at which the reflection should be predicted.
+ * Used by predict_to_res() and update_predictions()
+ **/
+void set_kpred(Reflection *refl, double kpred)
+{
+ refl->data.kpred = kpred;
+}
+
+
+/**
+ * set_exerr:
+ * @refl: A %Reflection
+ * @exerr: The excitation error for the reflection
*
**/
-void set_partial(Reflection *refl, double rlow, double rhigh, double p)
+void set_exerr(Reflection *refl, double exerr)
{
- refl->data.rlow = rlow;
- refl->data.rhigh = rhigh;
- refl->data.p = p;
+ refl->data.exerr = exerr;
}
diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h
index f2126ac9..9f494474 100644
--- a/libcrystfel/src/reflist.h
+++ b/libcrystfel/src/reflist.h
@@ -3,11 +3,11 @@
*
* Fast reflection/peak list
*
- * Copyright © 2012-2016 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2011-2016 Thomas White <taw@physics.org>
+ * 2011-2017 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -87,6 +87,9 @@ extern Reflection *next_found_refl(Reflection *refl);
extern void get_detector_pos(const Reflection *refl, double *fs, double *ss);
extern struct panel *get_panel(const Reflection *refl);
extern double get_partiality(const Reflection *refl);
+extern double get_khalf(const Reflection *refl);
+extern double get_kpred(const Reflection *refl);
+extern double get_exerr(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);
@@ -94,8 +97,6 @@ extern void get_symmetric_indices(const Reflection *refl,
signed int *hs, signed int *ks,
signed int *ls);
extern double get_intensity(const Reflection *refl);
-extern void get_partial(const Reflection *refl, double *rlow, double *rhigh,
- double *p);
extern int get_redundancy(const Reflection *refl);
extern double get_temp1(const Reflection *refl);
extern double get_temp2(const Reflection *refl);
@@ -109,7 +110,9 @@ extern int get_flag(const Reflection *refl);
extern void copy_data(Reflection *to, const Reflection *from);
extern void set_detector_pos(Reflection *refl, double fs, double ss);
extern void set_panel(Reflection *refl, struct panel *p);
-extern void set_partial(Reflection *refl, double rlow, double rhigh, double p);
+extern void set_kpred(Reflection *refl, double kpred);
+extern void set_khalf(Reflection *refl, double khalf);
+extern void set_exerr(Reflection *refl, double exerr);
extern void set_partiality(Reflection *refl, double p);
extern void set_lorentz(Reflection *refl, double L);
extern void set_intensity(Reflection *refl, double intensity);
diff --git a/scripts/plot-contourmap b/scripts/plot-contourmap
new file mode 100755
index 00000000..053a6be9
--- /dev/null
+++ b/scripts/plot-contourmap
@@ -0,0 +1,153 @@
+#!/usr/bin/env python
+
+import numpy as np
+import matplotlib
+import matplotlib.cm as cm
+import matplotlib.mlab as mlab
+import matplotlib.pyplot as plt
+import sys
+import fnmatch
+import os
+from matplotlib.widgets import RadioButtons,Button
+
+im=None
+cnt=None
+centre_marker=None
+
+def next_click(w):
+ global crystal
+ crystal += 20
+ print("Crystal %i" % crystal)
+ update_graph()
+
+
+def prev_click(w):
+ global crystal
+ crystal -= 20
+ print("Crystal %i" % crystal)
+ update_graph()
+
+
+def iteration_click(label):
+ global cycle
+ cycle = label.split(" ")[1].split("\n")[0]
+ update_graph()
+
+
+def variable_click(label):
+ global varpair
+ varpair = label
+ update_graph()
+
+
+def update_graph():
+
+ global im, cnt, centre_marker
+
+ filename="pr-logs/grid-crystal%i-cycle%s-%s.dat" % (crystal,cycle,varpair)
+ print filename
+
+ with open(filename, "r") as f:
+
+ l = f.readline().split(None,3)
+
+ min1 = float(l[0])
+ max1 = float(l[1])
+ cur1 = float(l[2])
+ label1 = l[3]
+
+ l = f.readline().split(None,3)
+ min2 = float(l[0])
+ max2 = float(l[1])
+ cur2 = float(l[2])
+ label2 = l[3]
+
+ extent = (min1, max1, min2, max2)
+ Z = np.loadtxt(fname=filename, ndmin=2, delimiter=" ", skiprows=2)
+
+ ax.set_xlim([min1,max1])
+ ax.set_ylim([min2,max2])
+
+ if not im:
+ im = ax.imshow(Z, interpolation='none', origin='lower',
+ cmap=cm.gray, extent=extent, aspect='auto')
+ im.autoscale()
+ else:
+ im.set_data(Z)
+ im.set_extent(extent)
+ im.autoscale()
+
+ levels = np.arange(0.0, 1.6, 0.1)
+
+ if cnt:
+ for coll in cnt.collections:
+ coll.remove()
+
+ cnt = ax.contour(Z, levels, origin='lower', linewidths=1, extent=extent)
+
+ ax.set_title(filename)
+ ax.set_xlabel(label1)
+ ax.set_ylabel(label2)
+ plt.flag()
+
+ if centre_marker:
+ centre_marker.remove()
+ centre_marker, = plt.plot(cur1, cur2, 'bH', color='r')
+
+ fig.canvas.draw()
+
+fig = plt.figure(figsize=(10,5))
+fig.subplots_adjust(left=0.05, bottom=0.05, right=0.70, top=0.95)
+
+# Find out what there is to plot
+crystals = []
+cycles = []
+varpairs = []
+for file in os.listdir("pr-logs"):
+ if not fnmatch.fnmatch(file, "grid-*.dat"):
+ continue
+ sp = file.rstrip(".dat").split("-")
+ crystals.append(int(sp[1].lstrip("crystal")))
+ cycles.append(sp[2].lstrip("cycle"))
+ varpairs.append(sp[3]+"-"+sp[4])
+
+crystals = sorted(set(crystals))
+cycles = sorted(set(cycles))
+varpairs = sorted(set(varpairs))
+
+crystal = crystals[0]
+cycle = cycles[0]
+varpair = varpairs[0]
+
+# Iteration selector
+ax = plt.axes([0.75, 0.55, 0.2, 0.40], facecolor="lightgoldenrodyellow")
+iterations = ["Iteration "+str(f) for f in cycles]
+if iterations[0] != "Iteration 0":
+ print("Whoops, couldn't find iteration 0!")
+ sys.exit(1)
+iterations[0] = "Iteration 0\n(initial scaling only)"
+if iterations[len(iterations)-1] != "Iteration F":
+ print("Whoops, couldn't find final iteration!")
+else:
+ iterations[len(iterations)-1] = "Iteration F\n(after final merge)"
+iteration = RadioButtons(ax, iterations)
+iteration.on_clicked(iteration_click)
+
+# Variable selector
+ax = plt.axes([0.75, 0.20, 0.2, 0.30], facecolor="lightgoldenrodyellow")
+variable = RadioButtons(ax, varpairs)
+variable.on_clicked(variable_click)
+
+# Crystal selector
+ax = plt.axes([0.75, 0.08, 0.20, 0.06])
+crystal_prev = Button(ax, "Previous crystal")
+crystal_prev.on_clicked(prev_click)
+ax = plt.axes([0.75, 0.01, 0.20, 0.06])
+crystal_next = Button(ax, "Next crystal")
+crystal_next.on_clicked(next_click)
+
+ax = plt.axes([0.1, 0.1, 0.6, 0.8])
+update_graph()
+plt.colorbar(im, shrink=0.8)
+
+plt.show()
diff --git a/scripts/plot-pr b/scripts/plot-pr
new file mode 100755
index 00000000..25427cae
--- /dev/null
+++ b/scripts/plot-pr
@@ -0,0 +1,196 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Plot post-refinement results
+#
+# Copyright © 2017-2018 Deutsches Elektronen-Synchrotron DESY,
+# a research centre of the Helmholtz Association.
+#
+# Author:
+# 2017-2018 Thomas White <taw@physics.org>
+#
+
+import sys
+import matplotlib
+import numpy as np
+from matplotlib import pyplot as plt
+from operator import truediv, itemgetter
+from matplotlib.widgets import RadioButtons,SpanSelector,Button
+
+crystal=0
+resmin=0
+resmax=100e9
+inum="0"
+first=True
+
+def update_graph():
+
+ global inum, resmin, resmax, first
+
+ # Update pobs/pcalc scatter plot
+ xlim = pgraphE.get_xlim()
+ ylim = pgraphE.get_ylim()
+ pgraphE.cla()
+ pgraphE.scatter([x[0] for x in pgraph if x[2]==inum if x[3]<resmax if x[3]>resmin],
+ [x[1] for x in pgraph if x[2]==inum if x[3]<resmax if x[3]>resmin], marker="+")
+ pgraphE.grid(b=True, which='major', color='k', linestyle='--')
+ if not first:
+ pgraphE.set_xlim(xlim)
+ pgraphE.set_ylim(ylim)
+ pgraphE.set_xlabel("pcalc")
+ pgraphE.set_ylabel("pobs")
+ pgraphE.set_title("Observed and calculated partialities")
+
+ # Update spectrum plot
+ pobsE.set_data([x[0] for x in specgraph if x[3]==inum if x[4]<resmax if x[4]>resmin],
+ [x[2] for x in specgraph if x[3]==inum if x[4]<resmax if x[4]>resmin])
+ pcalcE.set_data([x[0] for x in specgraph if x[3]==inum if x[4]<resmax if x[4]>resmin],
+ [x[1] for x in specgraph if x[3]==inum if x[4]<resmax if x[4]>resmin])
+ specgraphE.set_title("Spectrum plot for crystal %i" % crystal)
+
+ first = False
+ plt.draw()
+
+
+def iteration_click(label):
+ global inum
+ inum = label.split(" ")[1].split("\n")[0]
+ update_graph()
+
+
+def resolutionchange(rmin,rmax):
+ global resmin,resmax
+ resmin = rmin
+ resmax = rmax
+ update_graph()
+
+
+def read_spectrum(filename):
+
+ specgraph = []
+
+ try:
+ fh = open(filename, 'r')
+ except IOError:
+ print("Failed to read "+filename)
+ raise IOError
+
+ fline = fh.readline()
+ fline = fh.readline() # Ignore header (two lines)
+
+ while True:
+ fline = fh.readline()
+ if not fline:
+ break
+ kval = float(fline.split()[0])
+ res = float(fline.split()[1])
+ pcalc = float(fline.split()[2])
+ pobs = float(fline.split()[3])
+ inum = fline.split()[4]
+ specgraph.append([kval,pcalc,pobs,inum,res])
+
+ return sorted(specgraph, key=itemgetter(0))
+
+
+def next_click(w):
+ global specgraph, crystal
+ crystal += 20
+ specgraph = read_spectrum("pr-logs/specgraph-crystal%i.dat" % crystal)
+ update_graph()
+
+
+def prev_click(w):
+ global specgraph, crystal
+ crystal -= 20
+ specgraph = read_spectrum("pr-logs/specgraph-crystal%i.dat" % crystal)
+ update_graph()
+
+
+def set_sensible_axes(w):
+ pgraphE.set_xlim([-0.1,1.1])
+ pgraphE.set_ylim([-0.2,2.0])
+ specgraphE.set_ylim([-5.0,5.0])
+ plt.draw()
+
+
+# Read the pgraph file
+pgraph = []
+
+try:
+ fh = open("pr-logs/pgraph.dat", 'r')
+except IOError:
+ print("Failed to read "+filename)
+ raise IOError
+
+fline = fh.readline() # Ignore header
+
+while True:
+ fline = fh.readline()
+ if not fline:
+ break
+ res = float(fline.split()[4])
+ Ipart = float(fline.split()[5])
+ pc = float(fline.split()[6])
+ po = float(fline.split()[7])
+ inum = fline.split()[8]
+ pgraph.append([pc,po,inum,res,Ipart])
+
+
+# Read the spectrum file
+specgraph = read_spectrum("pr-logs/specgraph-crystal%i.dat" % crystal)
+
+fig = plt.figure(figsize=(15,6))
+fig.subplots_adjust(left=0.05, right=0.65)
+
+# pobs/pcalc
+pgraphE = fig.add_subplot(121)
+
+# Spectrum
+specgraphE = fig.add_subplot(122)
+pobsE = specgraphE.plot([0,1],[0,1],label="pobs")[0]
+pcalcE = specgraphE.plot([0,1],[0,1],label="pcalc")[0]
+specgraphE.set_xlabel("khalf / m^-1")
+specgraphE.set_ylabel("Partiality")
+specgraphE.grid(b=True, which='major', color='k', linestyle='--')
+specgraphE.legend()
+
+# Iteration selector
+ax = plt.axes([0.75, 0.50, 0.2, 0.40], facecolor="lightgoldenrodyellow")
+iterations = sorted(set(["Iteration "+str(x[2]) for x in pgraph]))
+if iterations[0] != "Iteration 0":
+ print("Whoops, couldn't find iteration 0!")
+ sys.exit(1)
+iterations[0] = "Iteration 0\n(initial scaling only)"
+if iterations[len(iterations)-1] != "Iteration F":
+ print("Whoops, couldn't find final iteration!")
+else:
+ iterations[len(iterations)-1] = "Iteration F\n(after final merge)"
+iteration = RadioButtons(ax, iterations)
+iteration.on_clicked(iteration_click)
+iteration.set_active(len(iterations)-1)
+
+# Resolution selector
+ax = plt.axes([0.75, 0.25, 0.2, 0.15], facecolor="lightgoldenrodyellow")
+ax.hist([x[3] for x in pgraph],weights=[x[4] for x in pgraph],bins=100)
+ax.set_xlabel("Resolution / m^-1")
+ax.set_ylabel("Frequency")
+ax.set_title("Resolution selector")
+resolution = SpanSelector(ax, resolutionchange, 'horizontal', span_stays=True, useblit=True)
+
+# Crystal switcher
+ax = plt.axes([0.75, 0.10, 0.09, 0.05])
+crystal_prev = Button(ax, "Previous crystal")
+crystal_prev.on_clicked(prev_click)
+ax = plt.axes([0.86, 0.10, 0.09, 0.05])
+crystal_next = Button(ax, "Next crystal")
+crystal_next.on_clicked(next_click)
+
+# Set sensible axes
+ax = plt.axes([0.75, 0.01, 0.09, 0.05])
+sensible_axes = Button(ax, "Set sensible axes")
+sensible_axes.on_clicked(set_sensible_axes)
+
+specgraphE.relim()
+specgraphE.autoscale_view(True,True,True)
+
+plt.show()
diff --git a/src/merge.c b/src/merge.c
index 17e57c90..e60aca15 100644
--- a/src/merge.c
+++ b/src/merge.c
@@ -50,7 +50,7 @@
/* Minimum partiality of a reflection for it to be merged */
-#define MIN_PART_MERGE (0.05)
+#define MIN_PART_MERGE (0.3)
struct merge_queue_args
@@ -59,7 +59,6 @@ struct merge_queue_args
pthread_rwlock_t full_lock;
Crystal **crystals;
int n_started;
- PartialityModel pmodel;
double push_res;
int use_weak;
long long int n_reflections;
@@ -159,7 +158,6 @@ static void run_merge_job(void *vwargs, int cookie)
signed int h, k, l;
double mean, sumweight, M2, temp, delta, R;
double corr, res, w;
- //double esd;
if ( get_partiality(refl) < MIN_PART_MERGE ) continue;
@@ -189,7 +187,7 @@ static void run_merge_job(void *vwargs, int cookie)
}
/* Total (multiplicative) correction factor */
- corr = exp(-G) * exp(B*res*res) * get_lorentz(refl)
+ corr = G * exp(B*res*res) * get_lorentz(refl)
/ get_partiality(refl);
if ( isnan(corr) ) {
ERROR("NaN corr:\n");
@@ -200,8 +198,8 @@ static void run_merge_job(void *vwargs, int cookie)
continue;
}
- //esd = get_esd_intensity(refl) * corr;
- w = 1.0;
+ /* Reflections count less the more they have to be scaled up */
+ w = 1.0 / corr;
/* Running mean and variance calculation */
temp = w + sumweight;
@@ -229,7 +227,7 @@ static void finalise_merge_job(void *vqargs, void *vwargs)
RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
- PartialityModel pmodel, int min_meas,
+ int min_meas,
double push_res, int use_weak)
{
RefList *full;
@@ -245,7 +243,6 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
qargs.full = full;
qargs.n_started = 0;
qargs.crystals = crystals;
- qargs.pmodel = pmodel;
qargs.push_res = push_res;
qargs.use_weak = use_weak;
qargs.n_reflections = 0;
diff --git a/src/merge.h b/src/merge.h
index 1911584f..56229e4c 100644
--- a/src/merge.h
+++ b/src/merge.h
@@ -40,7 +40,6 @@
#include "geometry.h"
extern RefList *merge_intensities(Crystal **crystals, int n, int n_threads,
- PartialityModel pmodel, int min_meas,
- double push_res, int use_weak);
+ int min_meas, double push_res, int use_weak);
#endif /* MERGE */
diff --git a/src/partial_sim.c b/src/partial_sim.c
index 4a348742..fc7845c1 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -386,7 +386,7 @@ static void run_job(void *vwargs, int cookie)
reflections = predict_to_res(cr, largest_q(&wargs->image));
crystal_set_reflections(cr, reflections);
- calculate_partialities(cr, PMODEL_SCSPHERE);
+ calculate_partialities(cr, PMODEL_XSPHERE);
for ( i=0; i<NBINS; i++ ) {
wargs->n_ref[i] = 0;
@@ -862,7 +862,7 @@ int main(int argc, char *argv[])
STATUS(" Background: none (no image "
"output)\n");
}
- STATUS(" Partiality model: scsphere (hardcoded)\n");
+ STATUS(" Partiality model: xsphere (hardcoded)\n");
STATUS(" Noise standard deviation: %.2f detector units\n",
noise_stddev);
if ( random_intensities ) {
@@ -874,6 +874,7 @@ int main(int argc, char *argv[])
input_file, sym_str);
}
STATUS(" Max error in cell components: %.2f %%\n", cnoise);
+ STATUS("Scale factor standard deviation: %.2f\n", osf_stddev);
if ( random_intensities ) {
full = reflist_new();
diff --git a/src/partialator.c b/src/partialator.c
index 428b0d68..21b02266 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -3,11 +3,11 @@
*
* Scaling and post refinement for coherent nanocrystallography
*
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2015 Thomas White <taw@physics.org>
+ * 2010-2018 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -40,6 +40,7 @@
#include <assert.h>
#include <pthread.h>
#include <gsl/gsl_errno.h>
+#include <sys/stat.h>
#include <image.h>
#include <utils.h>
@@ -190,7 +191,7 @@ static void write_split(Crystal **crystals, int n_crystals, const char *outfile,
}
snprintf(tmp, 1024, "%s1", outfile);
split = merge_intensities(crystals1, n_crystals1, nthreads,
- pmodel, min_measurements, push_res, 1);
+ min_measurements, push_res, 1);
if ( split == NULL ) {
ERROR("Not enough crystals for two way split!\n");
@@ -202,7 +203,7 @@ static void write_split(Crystal **crystals, int n_crystals, const char *outfile,
reflist_free(split);
snprintf(tmp, 1024, "%s2", outfile);
split = merge_intensities(crystals2, n_crystals2, nthreads,
- pmodel, min_measurements, push_res, 1);
+ min_measurements, push_res, 1);
STATUS("and %s\n", tmp);
write_reflist_2(tmp, split, sym);
reflist_free(split);
@@ -289,7 +290,7 @@ static void write_custom_split(struct custom_split *csplit, int dsn,
STATUS("Writing dataset '%s' to %s (%i crystals)\n",
csplit->dataset_names[dsn], tmp, n_crystalsn);
split = merge_intensities(crystalsn, n_crystalsn, nthreads,
- pmodel, min_measurements, push_res, 1);
+ min_measurements, push_res, 1);
write_reflist_2(tmp, split, sym);
reflist_free(split);
@@ -325,7 +326,12 @@ static void show_help(const char *s)
" -j <n> Run <n> analyses in parallel.\n"
" --no-free Disable cross-validation (testing only).\n"
" --custom-split List of files for custom dataset splitting.\n"
-" --max-rel-B Maximum allowable relative |B| factor.\n");
+" --max-rel-B Maximum allowable relative |B| factor.\n"
+" --no-logs Do not write extensive log files.\n"
+" -w <pg> Apparent point group for resolving ambiguities.\n"
+" --operator=<op> Indexing ambiguity operator for resolving.\n"
+" --force-bandwidth=<n> Set all bandwidths to <n> (fraction).\n"
+" --force-radius=<n> Set all profile radii to <n> nm^-1.\n");
}
@@ -579,7 +585,7 @@ static int set_initial_params(Crystal *cr, FILE *fh)
} else {
- crystal_set_osf(cr, 0.0);
+ crystal_set_osf(cr, 1.0);
crystal_set_Bfac(cr, 0.0);
}
@@ -604,13 +610,22 @@ static void select_free_reflections(RefList *list, gsl_rng *rng)
static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
- int fr)
+ int fr, signed int inum)
{
Reflection *refl;
RefListIterator *iter;
double G = crystal_get_osf(cr);
double B = crystal_get_Bfac(cr);
UnitCell *cell = crystal_get_cell(cr);
+ char ins[5];
+
+ if ( inum >= 0 ) {
+ snprintf(ins, 4, "%i", inum);
+ } else {
+ ins[0] = 'F';
+ ins[1] = '\0';
+ }
+
for ( refl = first_refl(list, &iter);
refl != NULL;
@@ -623,6 +638,9 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
if ( !get_flag(refl) ) continue; /* Not free-flagged */
+ /* Strong reflections only */
+ if ( get_intensity(refl) < 3.0*get_esd_intensity(refl) ) continue;
+
get_indices(refl, &h, &k, &l);
res = resolution(cell, h, k, l);
if ( 2.0*res > crystal_get_resolution_limit(cr) ) continue;
@@ -630,46 +648,55 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr,
match = find_refl(full, h, k, l);
if ( match == NULL ) continue;
+ /* Don't calculate pobs if reference reflection is weak */
+ if ( fabs(get_intensity(match)) / get_esd_intensity(match) < 3.0 ) continue;
+
/* Calculated partiality */
pcalc = get_partiality(refl);
/* Observed partiality */
- corr = exp(-G) * exp(B*res*res) * get_lorentz(refl);
+ corr = G * exp(B*res*res) * get_lorentz(refl);
Ipart = get_intensity(refl) * corr;
pobs = Ipart / get_intensity(match);
- fprintf(fh, "%5i %4i %4i %4i %8.4f %8.3f %8.3f\n",
- fr, h, k, l, 2*res/1e9, pcalc, pobs);
+ fprintf(fh, "%5i %4i %4i %4i %e %e %8.3f %8.3f %s\n",
+ fr, h, k, l, 2*res, Ipart, pcalc, pobs, ins);
}
}
static void write_pgraph(RefList *full, Crystal **crystals, int n_crystals,
- int iter)
+ signed int iter, const char *suff)
{
FILE *fh;
char tmp[256];
int i;
- snprintf(tmp, 256, "pgraph-iter%i.dat", iter);
+ snprintf(tmp, 256, "pr-logs/pgraph%s.dat", suff);
+
+ if ( iter == 0 ) {
+ fh = fopen(tmp, "w");
+ } else {
+ fh = fopen(tmp, "a");
+ }
- fh = fopen(tmp, "w");
if ( fh == NULL ) {
ERROR("Failed to open '%s'\n", tmp);
return;
}
- fprintf(fh, " fr h k l 1/d(nm) pcalc pobs\n");
+ if ( iter == 0 ) {
+ fprintf(fh, " Crystal h k l 1/d(m) Ipart pcalc pobs iteration\n");
+ }
for ( i=0; i<n_crystals; i++ ) {
if ( crystal_get_user_flag(crystals[i]) != 0 ) continue;
write_to_pgraph(fh, crystal_get_reflections(crystals[i]), full,
- crystals[i], i);
+ crystals[i], i, iter);
}
fclose(fh);
-
}
@@ -690,8 +717,8 @@ static void all_residuals(Crystal **crystals, int n_crystals, RefList *full,
if ( crystal_get_user_flag(crystals[i]) ) continue;
- r = residual(crystals[i], full, 0, NULL, NULL);
- free_r = residual(crystals[i], full, 1, NULL, NULL);
+ r = residual(crystals[i], full, 0, NULL, NULL, 1);
+ free_r = residual(crystals[i], full, 1, NULL, NULL, 1);
log_r = log_residual(crystals[i], full, 0, NULL, NULL);
free_log_r = log_residual(crystals[i], full, 1, NULL, NULL);
@@ -720,32 +747,77 @@ static void show_all_residuals(Crystal **crystals, int n_crystals,
}
-static void dump_parameters(const char *filename, Crystal **crystals,
- int n_crystals)
+struct log_qargs
{
- FILE *fh;
- fh = fopen(filename, "w");
- if ( fh == NULL ) {
- ERROR("Couldn't open partialator.params!\n");
- } else {
- fprintf(fh, " cr OSF relB div"
- " flag filename event\n");
- int i;
- for ( i=0; i<n_crystals; i++ ) {
- struct image *img;
- char *evt_str;
- img = crystal_get_image(crystals[i]);
- evt_str = get_event_string(img->event);
- fprintf(fh, "%4i %10.5f %10.2f %8.5e %-25s %s %s\n",
- i, crystal_get_osf(crystals[i]),
- crystal_get_Bfac(crystals[i])*1e20,
- crystal_get_image(crystals[i])->div,
- str_prflag(crystal_get_user_flag(crystals[i])),
- img->filename, evt_str);
- free(evt_str);
- }
- fclose(fh);
- }
+ int iter;
+ int next;
+ Crystal **crystals;
+ int n_crystals;
+ RefList *full;
+};
+
+
+struct log_args
+{
+ Crystal *cr;
+ RefList *full;
+ int iter;
+ int cnum;
+};
+
+
+static void *get_log_task(void *vp)
+{
+ struct log_qargs *qargs = vp;
+ struct log_args *task;
+
+ if ( qargs->next >= qargs->n_crystals ) return NULL;
+
+ task = malloc(sizeof(struct log_args));
+ if ( task == NULL ) return NULL;
+
+ task->cr = qargs->crystals[qargs->next];
+ task->full = qargs->full;
+ task->iter = qargs->iter;
+ task->cnum = qargs->next;
+
+ qargs->next += 20;
+ return task;
+}
+
+
+static void write_logs(void *vp, int cookie)
+{
+ struct log_args *args = vp;
+ write_specgraph(args->cr, args->full, args->iter, args->cnum);
+ write_gridscan(args->cr, args->full, args->iter, args->cnum);
+ write_test_logs(args->cr, args->full, args->iter, args->cnum);
+}
+
+
+static void done_log(void *qargs, void *vp)
+{
+ struct log_args *task = vp;
+ free(task);
+}
+
+
+static void write_logs_parallel(Crystal **crystals, int n_crystals,
+ RefList *full, int iter, int n_threads)
+{
+ struct log_qargs qargs;
+
+
+ qargs.iter = iter;
+ qargs.next = 0;
+ qargs.full = full;
+ qargs.crystals = crystals;
+ qargs.n_crystals = n_crystals;
+
+ STATUS("Writing logs...\n");
+ run_threads(n_threads, write_logs, get_log_task, done_log, &qargs,
+ n_crystals/20, 0, 0, 0);
+ STATUS("Done.\n");
}
@@ -756,9 +828,10 @@ int main(int argc, char *argv[])
char *outfile = NULL;
char *sym_str = NULL;
SymOpList *sym;
+ SymOpList *amb;
+ SymOpList *w_sym;
int nthreads = 1;
int i;
- struct image *images;
int n_iter = 10;
RefList *full;
int n_images = 0;
@@ -770,7 +843,7 @@ int main(int argc, char *argv[])
Stream *st;
Crystal **crystals;
char *pmodel_str = NULL;
- PartialityModel pmodel = PMODEL_SCSPHERE;
+ PartialityModel pmodel = PMODEL_XSPHERE;
int min_measurements = 2;
char *rval;
int polarisation = 1;
@@ -786,6 +859,13 @@ int main(int argc, char *argv[])
char *csplit_fn = NULL;
struct custom_split *csplit = NULL;
double max_B = 1e-18;
+ char *rfile = NULL;
+ RefList *reference = NULL;
+ int no_logs = 0;
+ char *w_sym_str = NULL;
+ char *operator = NULL;
+ double force_bandwidth = -1.0;
+ double force_radius = -1.0;
/* Long options */
const struct option longopts[] = {
@@ -806,8 +886,12 @@ int main(int argc, char *argv[])
{"push-res", 1, NULL, 5},
{"res-push", 1, NULL, 5}, /* compat */
{"custom-split", 1, NULL, 6},
- {"max-rel-B" , 1, NULL, 7},
- {"max-rel-b" , 1, NULL, 7}, /* compat */
+ {"max-rel-B", 1, NULL, 7},
+ {"max-rel-b", 1, NULL, 7}, /* compat */
+ {"reference", 1, NULL, 8}, /* ssshhh! */
+ {"operator", 1, NULL, 9},
+ {"force-bandwidth", 1, NULL, 10},
+ {"force-radius", 1, NULL, 11},
{"no-scale", 0, &no_scale, 1},
{"no-pr", 0, &no_pr, 1},
@@ -817,6 +901,7 @@ int main(int argc, char *argv[])
{"polarization", 0, &polarisation, 1}, /* compat */
{"no-free", 0, &no_free, 1},
{"output-every-cycle", 0, &output_everycycle, 1},
+ {"no-logs", 0, &no_logs, 1},
{0, 0, NULL, 0}
};
@@ -828,7 +913,7 @@ int main(int argc, char *argv[])
}
/* Short options */
- while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:j:m:",
+ while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:j:m:w:",
longopts, NULL)) != -1)
{
@@ -887,6 +972,10 @@ int main(int argc, char *argv[])
pmodel_str = strdup(optarg);
break;
+ case 'w' :
+ w_sym_str = strdup(optarg);
+ break;
+
case 2 :
errno = 0;
min_measurements = strtod(optarg, &rval);
@@ -933,6 +1022,33 @@ int main(int argc, char *argv[])
max_B = max_B * 1e-20;
break;
+ case 8 :
+ rfile = strdup(optarg);
+ break;
+
+ case 9 :
+ operator = strdup(optarg);
+ break;
+
+ case 10 :
+ errno = 0;
+ force_bandwidth = strtod(optarg, &rval);
+ if ( *rval != '\0' ) {
+ ERROR("Invalid value for --force-bandwidth.\n");
+ return 1;
+ }
+ break;
+
+ case 11 :
+ errno = 0;
+ force_radius = strtod(optarg, &rval);
+ if ( *rval != '\0' ) {
+ ERROR("Invalid value for --force-radius.\n");
+ return 1;
+ }
+ force_radius *= 1e9;
+ break;
+
case 0 :
break;
@@ -972,13 +1088,49 @@ int main(int argc, char *argv[])
free(sym_str);
if ( sym == NULL ) return 1;
+ if ( (w_sym_str != NULL) && (operator != NULL) ) {
+ ERROR("Specify the apparent symmetry (-w) or the operator, "
+ "not both.\n");
+ return 1;
+ }
+
+ if ( w_sym_str == NULL ) {
+ w_sym = NULL;
+ amb = NULL;
+ } else {
+ pointgroup_warning(w_sym_str);
+ w_sym = get_pointgroup(w_sym_str);
+ free(w_sym_str);
+ if ( w_sym == NULL ) return 1;
+ amb = get_ambiguities(w_sym, sym);
+ if ( amb == NULL ) {
+ ERROR("Couldn't find ambiguity operator.\n");
+ ERROR("Check that your values for -y and -w are "
+ "correct.\n");
+ return 1;
+ }
+
+ }
+
+ if ( operator ) {
+ amb = parse_symmetry_operations(operator);
+ if ( amb == NULL ) return 1;
+ set_symmetry_name(amb, "Ambiguity");
+ }
+
+ if ( amb != NULL ) {
+ STATUS("Indexing ambiguity resolution enabled. "
+ "The ambiguity operation(s) are:\n");
+ describe_symmetry(amb);
+ /* In contrast to ambigator, partialator can deal with multiple
+ * ambiguities at once */
+ }
+
if ( pmodel_str != NULL ) {
if ( strcmp(pmodel_str, "unity") == 0 ) {
pmodel = PMODEL_UNITY;
- } else if ( strcmp(pmodel_str, "scgaussian") == 0 ) {
- pmodel = PMODEL_SCGAUSSIAN;
- } else if ( strcmp(pmodel_str, "scsphere") == 0 ) {
- pmodel = PMODEL_SCSPHERE;
+ } else if ( strcmp(pmodel_str, "xsphere") == 0 ) {
+ pmodel = PMODEL_XSPHERE;
} else if ( strcmp(pmodel_str, "random") == 0 ) {
pmodel = PMODEL_RANDOM;
} else {
@@ -989,7 +1141,22 @@ int main(int argc, char *argv[])
if ( (pmodel == PMODEL_UNITY) && !no_pr ) {
no_pr = 1;
- STATUS("Setting --no-pr because we are not modelling partialities (--model=unity).\n");
+ STATUS("Setting --no-pr because we are not modelling "
+ "partialities (--model=unity).\n");
+ }
+
+ if ( !no_logs ) {
+ int r = mkdir("pr-logs", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+ if ( r ) {
+ if ( errno == EEXIST ) {
+ ERROR("A folder called 'pr-logs' exists in the "
+ "working directory.\n");
+ ERROR("Please delete or move it first.\n");
+ } else {
+ ERROR("Failed to create pr-logs folder.\n");
+ }
+ return 1;
+ }
}
/* Read the custom split list (if applicable) */
@@ -1002,6 +1169,20 @@ int main(int argc, char *argv[])
free(csplit_fn);
}
+ if ( rfile != NULL ) {
+ RefList *rread;
+ rread = read_reflections(rfile);
+ if ( rread == NULL ) {
+ ERROR("Failed to read reference reflections\n");
+ return 1;
+ }
+ reference = asymmetric_indices(rread, sym);
+ reflist_free(rread);
+ ERROR("WARNING: Using an external reference.\n");
+ ERROR("WARNING: If you publish a structure based on the result,"
+ " expect to have to retract your paper!\n");
+ }
+
gsl_set_error_handler_off();
rng = gsl_rng_alloc(gsl_rng_mt19937);
@@ -1009,7 +1190,6 @@ int main(int argc, char *argv[])
n_images = 0;
n_crystals = 0;
n_crystals_seen = 0;
- images = NULL;
crystals = NULL;
if ( sparams_fn != NULL ) {
char line[1024];
@@ -1033,35 +1213,27 @@ int main(int argc, char *argv[])
RefList *as;
int i;
- struct image *images_new;
- struct image *cur;
+ struct image cur;
- images_new = realloc(images, (n_images+1)*sizeof(struct image));
- if ( images_new == NULL ) {
- ERROR("Failed to allocate memory for image list.\n");
- return 1;
- }
- images = images_new;
- cur = &images[n_images];
-
- cur->div = NAN;
- cur->bw = NAN;
- cur->det = NULL;
- if ( read_chunk_2(st, cur, STREAM_READ_REFLECTIONS
- | STREAM_READ_UNITCELL) != 0 ) {
+ cur.div = NAN;
+ cur.bw = NAN;
+ cur.det = NULL;
+ if ( read_chunk_2(st, &cur, STREAM_READ_REFLECTIONS
+ | STREAM_READ_UNITCELL) != 0 ) {
break;
}
- if ( isnan(cur->div) || isnan(cur->bw) ) {
+ if ( isnan(cur.div) || isnan(cur.bw) ) {
ERROR("Chunk doesn't contain beam parameters.\n");
return 1;
}
- for ( i=0; i<cur->n_crystals; i++ ) {
+ for ( i=0; i<cur.n_crystals; i++ ) {
Crystal *cr;
Crystal **crystals_new;
RefList *cr_refl;
+ struct image *image;
n_crystals_seen++;
if ( n_crystals_seen <= start_after ) continue;
@@ -1074,11 +1246,19 @@ int main(int argc, char *argv[])
return 1;
}
crystals = crystals_new;
- crystals[n_crystals] = cur->crystals[i];
+ crystals[n_crystals] = cur.crystals[i];
cr = crystals[n_crystals];
- /* Image pointer will change due to later reallocs */
- crystal_set_image(cr, NULL);
+ image = malloc(sizeof(struct image));
+ if ( image == NULL ) {
+ ERROR("Failed to allocatea memory for image.\n");
+ return 1;
+ }
+
+ crystal_set_image(cr, image);
+ *image = cur;
+ image->n_crystals = 1;
+ image->crystals = &crystals[n_crystals];
/* This is the raw list of reflections */
cr_refl = crystal_get_reflections(cr);
@@ -1088,7 +1268,7 @@ int main(int argc, char *argv[])
if ( polarisation ) {
polarisation_correction(cr_refl,
crystal_get_cell(cr),
- cur);
+ image);
}
if ( !no_free ) select_free_reflections(cr_refl, rng);
@@ -1124,65 +1304,76 @@ int main(int argc, char *argv[])
close_stream(st);
- /* Fill in image pointers */
- for ( i=0; i<n_images; i++ ) {
- int j;
- for ( j=0; j<images[i].n_crystals; j++ ) {
-
- Crystal *cryst;
-
- cryst = images[i].crystals[j];
- crystal_set_image(cryst, &images[i]);
-
- /* Now it's safe to do the following */
- update_predictions(cryst);
- calculate_partialities(cryst, pmodel);
-
+ STATUS("Initial partiality calculation...\n");
+ for ( i=0; i<n_crystals; i++ ) {
+ Crystal *cr = crystals[i];
+ if ( force_bandwidth > 0.0 ) {
+ crystal_get_image(cr)->bw = force_bandwidth;
+ }
+ if ( force_radius > 0.0 ) {
+ crystal_set_profile_radius(cr, force_radius);
}
+ update_predictions(cr);
+ calculate_partialities(cr, pmodel);
}
if (csplit != NULL) check_csplit(crystals, n_crystals, csplit);
/* Make a first pass at cutting out crap */
- STATUS("Checking patterns.\n");
+ STATUS("Early rejection...\n");
//early_rejection(crystals, n_crystals);
- /* Initial rejection, figures of merit etc */
- full = merge_intensities(crystals, n_crystals, nthreads, pmodel,
- min_measurements, push_res, 1);
+ /* Create reference data set if we don't already have one */
+ if ( reference == NULL ) {
+ STATUS("Initial scaling...\n");
+ scale_all(crystals, n_crystals, nthreads);
+ full = merge_intensities(crystals, n_crystals, nthreads,
+ min_measurements, push_res, 1);
+ } else {
+ full = reference;
+ }
+
+ /* Scale everything to the reference */
+ STATUS("Initial scaling to reference...\n");
+ if ( !no_scale ) {
+ scale_all_to_reference(crystals, n_crystals, full, nthreads);
+ }
+
+ /* Check rejection and write figures of merit */
check_rejection(crystals, n_crystals, full, max_B);
show_all_residuals(crystals, n_crystals, full);
- write_pgraph(full, crystals, n_crystals, 0);
+ if ( !no_pr && !no_logs ) {
+ write_pgraph(full, crystals, n_crystals, 0, "");
+ write_logs_parallel(crystals, n_crystals, full, 0, nthreads);
+ }
/* Iterate */
for ( i=0; i<n_iter; i++ ) {
STATUS("Scaling and refinement cycle %i of %i\n", i+1, n_iter);
- if ( !no_scale ) {
- scale_all(crystals, n_crystals, nthreads, pmodel);
- reflist_free(full);
- full = merge_intensities(crystals, n_crystals, nthreads,
- pmodel, min_measurements,
- push_res, 1);
+ if ( !no_pr ) {
+ refine_all(crystals, n_crystals, full, nthreads, pmodel,
+ 0, i+1, no_logs, sym, amb);
+ } else if ( !no_scale ) {
+ scale_all_to_reference(crystals, n_crystals, full, nthreads);
}
- if ( !no_pr ) {
- refine_all(crystals, n_crystals, full, nthreads,
- pmodel);
+ /* Create new reference if needed */
+ if ( reference == NULL ) {
reflist_free(full);
+ scale_all(crystals, n_crystals, nthreads);
full = merge_intensities(crystals, n_crystals, nthreads,
- pmodel, min_measurements,
+ min_measurements,
push_res, 1);
- }
+ } /* else full still equals reference */
check_rejection(crystals, n_crystals, full, max_B);
- reflist_free(full);
- full = merge_intensities(crystals, n_crystals, nthreads,
- pmodel, min_measurements,
- push_res, 1);
show_all_residuals(crystals, n_crystals, full);
- write_pgraph(full, crystals, n_crystals, i+1);
+
+ if ( !no_logs ) {
+ write_pgraph(full, crystals, n_crystals, i+1, "");
+ }
if ( output_everycycle ) {
@@ -1209,17 +1400,29 @@ int main(int argc, char *argv[])
}
}
- /* Dump parameters */
- snprintf(tmp, 1024, "iter%.2d_partialator.params", i+1);
- dump_parameters(tmp, crystals, n_crystals);
-
}
}
- full = merge_intensities(crystals, n_crystals, nthreads, pmodel,
- min_measurements, push_res, 1);
+ /* Final merge */
+ STATUS("Final merge...\n");
+ if ( reference == NULL ) {
+ reflist_free(full);
+ scale_all(crystals, n_crystals, nthreads);
+ full = merge_intensities(crystals, n_crystals, nthreads,
+ min_measurements,
+ push_res, 1);
+ } else {
+ scale_all_to_reference(crystals, n_crystals, reference, nthreads);
+ full = merge_intensities(crystals, n_crystals, nthreads,
+ min_measurements, push_res, 1);
+ }
+
+ /* Write final figures of merit (no rejection any more) */
show_all_residuals(crystals, n_crystals, full);
- write_pgraph(full, crystals, n_crystals, n_iter+1);
+ if ( !no_pr && !no_logs ) {
+ write_pgraph(full, crystals, n_crystals, -1, "");
+ write_logs_parallel(crystals, n_crystals, full, -1, nthreads);
+ }
/* Output results */
STATUS("Writing overall results to %s\n", outfile);
@@ -1239,9 +1442,6 @@ int main(int argc, char *argv[])
}
}
- /* Dump parameters */
- dump_parameters("partialator.params", crystals, n_crystals);
-
/* Clean up */
gsl_rng_free(rng);
for ( i=0; i<n_crystals; i++ ) {
@@ -1249,13 +1449,9 @@ int main(int argc, char *argv[])
crystal_free(crystals[i]);
}
reflist_free(full);
- free(sym);
+ free_symoplist(sym);
free(outfile);
free(crystals);
- for ( i=0; i<n_images; i++ ) {
- free(images[i].filename);
- }
- free(images);
free(infile);
return 0;
diff --git a/src/post-refinement.c b/src/post-refinement.c
index bf30d299..79ea187a 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -3,11 +3,11 @@
*
* Post refinement
*
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2015 Thomas White <taw@physics.org>
+ * 2010-2018 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -33,11 +33,8 @@
#include <stdlib.h>
#include <assert.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_vector.h>
-#include <gsl/gsl_linalg.h>
-#include <gsl/gsl_eigen.h>
-#include <gsl/gsl_blas.h>
+#include <gsl/gsl_multimin.h>
+#include <gsl/gsl_fit.h>
#include "image.h"
#include "post-refinement.h"
@@ -46,15 +43,13 @@
#include "geometry.h"
#include "cell.h"
#include "cell-utils.h"
+#include "reflist-utils.h"
+#include "scaling.h"
-/* Maximum number of iterations of NLSq to do for each image per macrocycle. */
-#define MAX_CYCLES (10)
-
struct prdata
{
int refined;
- int n_filtered;
};
const char *str_prflag(enum prflag flag)
@@ -79,625 +74,867 @@ const char *str_prflag(enum prflag flag)
case PRFLAG_BIGB :
return "B too big";
+ case PRFLAG_SCALEBAD :
+ return "bad scaling";
+
default :
return "Unknown flag";
}
}
-/* Returns dp(gauss)/dr at "r" */
-static double gaussian_fraction_gradient(double r, double R)
+double residual(Crystal *cr, const RefList *full, int free,
+ int *pn_used, const char *filename, int complain)
{
- const double ng = 2.6;
- const double sig = R/ng;
+ Reflection *refl;
+ RefListIterator *iter;
+ int n_used = 0;
+ double num = 0.0;
+ double den = 0.0;
+ double G = crystal_get_osf(cr);
+ double B = crystal_get_Bfac(cr);
+ UnitCell *cell = crystal_get_cell(cr);
- /* If the Ewald sphere isn't within the profile, the gradient is zero */
- if ( r < -R ) return 0.0;
- if ( r > +R ) return 0.0;
+ if ( linear_scale(full, crystal_get_reflections(cr), &G, complain) ) {
+ return GSL_NAN;
+ }
- return exp(-pow(r/sig, 2.0)/2.0) / (sig*sqrt(2.0*M_PI));
-}
+ for ( refl = first_refl(crystal_get_reflections(cr), &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ double p, w, corr, res;
+ signed int h, k, l;
+ Reflection *match;
+ double I_full;
+ double int1, int2;
+ if ( free && !get_flag(refl) ) continue;
-/* Returns dp(sph)/dr at "r" */
-static double sphere_fraction_gradient(double r, double pr)
-{
- double q, dpdq, dqdr;
+ get_indices(refl, &h, &k, &l);
+ res = resolution(cell, h, k, l);
+ match = find_refl(full, h, k, l);
+ if ( match == NULL ) continue;
+ I_full = get_intensity(match);
- /* If the Ewald sphere isn't within the profile, the gradient is zero */
- if ( r < -pr ) return 0.0;
- if ( r > +pr ) return 0.0;
+ if ( get_redundancy(match) < 2 ) continue;
- q = (r + pr)/(2.0*pr);
- dpdq = 6.0*(q - q*q);
- dqdr = 1.0 / (2.0*pr);
- return dpdq * dqdr;
-}
+ p = get_partiality(refl);
+ //if ( p < 0.2 ) continue;
+ corr = G * exp(B*res*res) * get_lorentz(refl);
+ int1 = get_intensity(refl) * corr;
+ int2 = p*I_full;
+ w = p;
-/* Returns dp/dr at "r" */
-static double partiality_gradient(double r, double pr,
- PartialityModel pmodel,
- double rlow, double rhigh)
-{
- double A, D;
+ num += fabs(int1 - int2) * w;
+ den += fabs(int1 + int2) * w/2.0;
- D = rlow - rhigh;
+ n_used++;
- switch ( pmodel ) {
+ }
+
+ if ( pn_used != NULL ) *pn_used = n_used;
+ return num/(den*sqrt(2));
+}
- default:
- case PMODEL_UNITY:
- return 0.0;
- case PMODEL_SCSPHERE:
- A = sphere_fraction_gradient(r, pr)/D;
- return 4.0*pr*A/3.0;
+static UnitCell *rotate_cell_xy(const UnitCell *cell, double ang1, double ang2)
+{
+ UnitCell *o;
+ double asx, asy, asz;
+ double bsx, bsy, bsz;
+ double csx, csy, csz;
+ double xnew, ynew, znew;
+
+ o = cell_new_from_cell(cell);
+
+ cell_get_reciprocal(o, &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+
+ /* "a" around x */
+ xnew = asx;
+ ynew = asy*cos(ang1) + asz*sin(ang1);
+ znew = -asy*sin(ang1) + asz*cos(ang1);
+ asx = xnew; asy = ynew; asz = znew;
+
+ /* "b" around x */
+ xnew = bsx;
+ ynew = bsy*cos(ang1) + bsz*sin(ang1);
+ znew = -bsy*sin(ang1) + bsz*cos(ang1);
+ bsx = xnew; bsy = ynew; bsz = znew;
+
+ /* "c" around x */
+ xnew = csx;
+ ynew = csy*cos(ang1) + csz*sin(ang1);
+ znew = -csy*sin(ang1) + csz*cos(ang1);
+ csx = xnew; csy = ynew; csz = znew;
+
+ /* "a" around y */
+ xnew = asx*cos(ang2) + asz*sin(ang2);
+ ynew = asy;
+ znew = -asx*sin(ang2) + asz*cos(ang2);
+ asx = xnew; asy = ynew; asz = znew;
+
+ /* "b" around y */
+ xnew = bsx*cos(ang2) + bsz*sin(ang2);
+ ynew = bsy;
+ znew = -bsx*sin(ang2) + bsz*cos(ang2);
+ bsx = xnew; bsy = ynew; bsz = znew;
+
+ /* "c" around y */
+ xnew = csx*cos(ang2) + csz*sin(ang2);
+ ynew = csy;
+ znew = -csx*sin(ang2) + csz*cos(ang2);
+ csx = xnew; csy = ynew; csz = znew;
+
+ cell_set_reciprocal(o, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
+
+ return o;
+}
- case PMODEL_SCGAUSSIAN:
- A = gaussian_fraction_gradient(r, pr)/D;
- return 4.0*pr*A/3.0;
+static const char *get_label(enum gparam p)
+{
+ switch ( p ) {
+ case GPARAM_ANG1 : return "First angle (radians)";
+ case GPARAM_ANG2 : return "Second angle (radians)";
+ case GPARAM_R : return "Profile radius (m^-1)";
+ case GPARAM_WAVELENGTH : return "Wavelength (Angstroms)";
+ default : return "unknown";
}
}
-static double sphere_fraction_rgradient(double r, double R)
+/* We set all the step sizes to 1, then scale them.
+ * This way, the size of the simplex stays meaningful and we possibly also
+ * avoid some roundoff errors */
+static double get_scale(enum gparam p)
{
- /* If the Ewald sphere isn't within the profile, the gradient is zero */
- if ( r < -R ) return 0.0;
- if ( r > +R ) return 0.0;
+ switch ( p ) {
+
+ case GPARAM_ANG1 : return deg2rad(0.05);
+ case GPARAM_ANG2 : return deg2rad(0.05);
+ case GPARAM_R : return 0.0005e9;
+ case GPARAM_WAVELENGTH : return 0.001e-10;
- return -(3.0*r/(4.0*R*R)) * (1.0 - r*r/(R*R));
+ default : return 0.0;
+
+ }
}
-static double gaussian_fraction_rgradient(double r, double R)
-{
- const double ng = 2.6;
- const double sig = R/ng;
+struct rf_priv {
+ const Crystal *cr;
+ const RefList *full;
+ enum gparam rv[32];
+ int verbose;
+ int serial;
+ gsl_vector *initial;
- /* If the Ewald sphere isn't within the profile, the gradient is zero */
- if ( r < -R ) return 0.0;
- if ( r > +R ) return 0.0;
+ /* For freeing later */
+ gsl_vector *vals;
+ gsl_vector *step;
- return -(ng*r/(sqrt(2.0*M_PI)*R*R))*exp(-r*r/(2.0*sig*sig));
+ /* So that it stays around until the end of minimisation */
+ gsl_multimin_function f;
+};
+
+
+static double get_actual_val(const gsl_vector *v, const gsl_vector *initial,
+ enum gparam *rv, int i)
+{
+ return gsl_vector_get(v, i) * get_scale(rv[i])
+ + gsl_vector_get(initial, i);
}
-static double volume_fraction_rgradient(double r, double pr,
- PartialityModel pmodel)
+static void apply_parameters(const gsl_vector *v, const gsl_vector *initial,
+ enum gparam *rv, Crystal *cr)
{
- switch ( pmodel )
- {
- case PMODEL_UNITY :
- return 1.0;
+ int i;
+ double ang1, ang2, R, lambda;
+ UnitCell *cell;
- case PMODEL_SCSPHERE :
- return sphere_fraction_rgradient(r, pr);
+ /* Default parameters if not used in refinement */
+ ang1 = 0.0;
+ ang2 = 0.0;
+ R = crystal_get_profile_radius(cr);
+ lambda = crystal_get_image(cr)->lambda;
- case PMODEL_SCGAUSSIAN :
- return gaussian_fraction_rgradient(r, pr);
+ for ( i=0; i<v->size; i++ ) {
- default :
- ERROR("No pmodel in volume_fraction_rgradient!\n");
- return 1.0;
- }
-}
+ double val;
+ val = get_actual_val(v, initial, rv, i);
-static double volume_fraction(double rlow, double rhigh, double pr,
- PartialityModel pmodel)
-{
- switch ( pmodel )
- {
- case PMODEL_UNITY :
- return 1.0;
+ switch ( rv[i] ) {
- case PMODEL_SCSPHERE :
- return sphere_fraction(rlow, rhigh, pr);
+ case GPARAM_ANG1 :
+ ang1 = val;
+ break;
- case PMODEL_SCGAUSSIAN :
- return gaussian_fraction(rlow, rhigh, pr);
+ case GPARAM_ANG2 :
+ ang2 = val;
+ break;
- default :
- ERROR("No pmodel in volume_fraction!\n");
- return 1.0;
- }
-}
+ case GPARAM_R :
+ R = val;
+ break;
+ case GPARAM_WAVELENGTH :
+ lambda = val;
+ break;
-/* Return the gradient of "fx" wrt parameter 'k' given the current
- * status of the crystal. */
-double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel)
-{
- double glow, ghigh;
- double rlow, rhigh, p;
- struct image *image = crystal_get_image(cr);
- double R = crystal_get_profile_radius(cr);
- double gr;
+ default :
+ ERROR("Don't understand parameter %i\n", rv[i]);
+ break;
- get_partial(refl, &rlow, &rhigh, &p);
+ }
+ }
- if ( k == GPARAM_R ) {
+ cell = rotate_cell_xy(crystal_get_cell_const(cr), ang1, ang2);
+ crystal_set_cell(cr, cell);
- double Rglow, Rghigh;
- double D, psph;
+ crystal_set_profile_radius(cr, R);
+ crystal_get_image(cr)->lambda = lambda;
+}
- D = rlow - rhigh;
- psph = volume_fraction(rlow, rhigh, R, pmodel);
- Rglow = volume_fraction_rgradient(rlow, R, pmodel);
- Rghigh = volume_fraction_rgradient(rhigh, R, pmodel);
+static double residual_f(const gsl_vector *v, void *pp)
+{
+ struct rf_priv *pv = pp;
+ RefList *list;
+ struct image im;
+ Crystal *cr;
+ UnitCell *cell;
+ double res;
+ int i;
- gr = 4.0*psph/(3.0*D) + (4.0*R/(3.0*D))*(Rglow - Rghigh);
- return gr;
+ for ( i=0; i<v->size; i++ ) {
+ if ( gsl_vector_get(v, i) > 100.0 ) return GSL_NAN;
+ }
+ cr = crystal_copy(pv->cr);
+ cell = cell_new_from_cell(crystal_get_cell(cr));
+ if ( cell == NULL ) return GSL_NAN;
+ crystal_set_cell(cr, cell);
+ im = *crystal_get_image(cr);
+ crystal_set_image(cr, &im);
+ apply_parameters(v, pv->initial, pv->rv, cr);
+
+ if ( fabs(crystal_get_profile_radius(cr)) > 5e9 ) {
+ crystal_free(cr);
+ if ( pv->verbose ) STATUS("radius > 5e9\n");
+ return GSL_NAN;
}
- /* Calculate the gradient of partiality wrt excitation error. */
- glow = partiality_gradient(rlow, R, pmodel, rlow, rhigh);
- ghigh = partiality_gradient(rhigh, R, pmodel, rlow, rhigh);
+ /* Can happen with grid scans and certain --force-radius values */
+ if ( fabs(crystal_get_profile_radius(cr)) < 0.0000001e9 ) {
+ crystal_free(cr);
+ if ( pv->verbose ) STATUS("radius very small\n");
+ return GSL_NAN;
+ }
- if ( k == GPARAM_DIV ) {
+ if ( im.lambda <= 0.0 ) {
+ crystal_free(cr);
+ if ( pv->verbose ) STATUS("lambda < 0\n");
+ return GSL_NAN;
+ }
- double D, psph, ds;
- signed int hs, ks, ls;
+ list = copy_reflist(crystal_get_reflections(cr));
+ crystal_set_reflections(cr, list);
- D = rlow - rhigh;
- psph = volume_fraction(rlow, rhigh, R, pmodel);
- get_symmetric_indices(refl, &hs, &ks, &ls);
- ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls);
+ update_predictions(cr);
+ calculate_partialities(cr, PMODEL_XSPHERE);
- gr = (ds/2.0)*(glow+ghigh) - 4.0*R*psph*ds/(3.0*D*D);
- return gr;
+ res = residual(cr, pv->full, 0, NULL, NULL, 0);
- }
+ cell_free(crystal_get_cell(cr));
+ reflist_free(crystal_get_reflections(cr));
+ crystal_free(cr);
- gr = r_gradient(crystal_get_cell(cr), k, refl, image) * (glow-ghigh);
- return gr;
+ return res;
}
-static void apply_cell_shift(UnitCell *cell, int k, double shift)
+static double get_initial_param(Crystal *cr, enum gparam p)
{
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
- double as, bs, cs;
+ switch ( p ) {
- cell_get_reciprocal(cell, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz);
- as = modulus(asx, asy, asz);
- bs = modulus(bsx, bsy, bsz);
- cs = modulus(csx, csy, csz);
+ case GPARAM_ANG1 : return 0.0;
+ case GPARAM_ANG2 : return 0.0;
+ case GPARAM_R : return crystal_get_profile_radius(cr);
+ case GPARAM_WAVELENGTH : return crystal_get_image(cr)->lambda;
- /* Forbid any step which looks too large */
- switch ( k )
- {
- case GPARAM_ASX :
- case GPARAM_ASY :
- case GPARAM_ASZ :
- if ( fabs(shift) > 0.1*as ) return;
- break;
-
- case GPARAM_BSX :
- case GPARAM_BSY :
- case GPARAM_BSZ :
- if ( fabs(shift) > 0.1*bs ) return;
- break;
-
- case GPARAM_CSX :
- case GPARAM_CSY :
- case GPARAM_CSZ :
- if ( fabs(shift) > 0.1*cs ) return;
- break;
- }
-
- switch ( k )
- {
- case GPARAM_ASX : asx += shift; break;
- case GPARAM_ASY : asy += shift; break;
- case GPARAM_ASZ : asz += shift; break;
- case GPARAM_BSX : bsx += shift; break;
- case GPARAM_BSY : bsy += shift; break;
- case GPARAM_BSZ : bsz += shift; break;
- case GPARAM_CSX : csx += shift; break;
- case GPARAM_CSY : csy += shift; break;
- case GPARAM_CSZ : csz += shift; break;
- }
-
- cell_set_reciprocal(cell, asx, asy, asz,
- bsx, bsy, bsz,
- csx, csy, csz);
+ default: return 0.0;
+
+ }
}
-/* Apply the given shift to the 'k'th parameter of 'image'. */
-static void apply_shift(Crystal *cr, int k, double shift)
+static int check_angle_shifts(gsl_vector *cur, const gsl_vector *initial,
+ enum gparam *rv, struct rf_priv *residual_f_priv)
{
- double t;
- struct image *image = crystal_get_image(cr);
-
- switch ( k ) {
-
- case GPARAM_DIV :
- if ( shift > 0.1*image->div ) return;
- image->div += shift;
- if ( image->div < 0.0 ) image->div = 0.0;
- break;
-
- case GPARAM_R :
- t = crystal_get_profile_radius(cr);
- if ( shift > 0.1*t ) return;
- t += shift;
- crystal_set_profile_radius(cr, t);
- break;
-
- case GPARAM_ASX :
- case GPARAM_ASY :
- case GPARAM_ASZ :
- case GPARAM_BSX :
- case GPARAM_BSY :
- case GPARAM_BSZ :
- case GPARAM_CSX :
- case GPARAM_CSY :
- case GPARAM_CSZ :
- apply_cell_shift(crystal_get_cell(cr), k, shift);
- break;
+ int i = 0;
+ double ang = 0.0;
- default :
- ERROR("No shift defined for parameter %i\n", k);
- abort();
+ while ( rv[i] != GPARAM_EOL ) {
+ if ( (rv[i] == GPARAM_ANG1) || (rv[i] == GPARAM_ANG2) ) {
+ ang += fabs(get_actual_val(cur, initial, rv, i));
+ }
+ rv++;
+ }
+ if ( rad2deg(ang) > 5.0 ) {
+ ERROR("More than 5 degrees total rotation!\n");
+ residual_f_priv->verbose = 1;
+ double res = residual_f(cur, residual_f_priv);
+ STATUS("residual after rotation = %e\n", res);
+ residual_f_priv->verbose = 2;
+ res = residual_f(initial, residual_f_priv);
+ STATUS("residual before rotation = %e\n", res);
+ return 1;
}
+ return 0;
}
-/* Perform one cycle of post refinement on 'image' against 'full' */
-static double pr_iterate(Crystal *cr, const RefList *full,
- PartialityModel pmodel,
- int *n_filtered, int verbose)
+static RefList *reindex_reflections(RefList *input, SymOpList *amb,
+ SymOpList *sym, int idx)
{
- gsl_matrix *M;
- gsl_vector *v;
- gsl_vector *shifts;
- int param;
+ RefList *n;
Reflection *refl;
RefListIterator *iter;
- RefList *reflections;
- double max_shift;
- int nref = 0;
- int num_params = 0;
- enum gparam rv[32];
- double G, B;
-
- if ( n_filtered != NULL ) *n_filtered = 0;
-
- rv[num_params++] = GPARAM_ASX;
- rv[num_params++] = GPARAM_ASY;
- rv[num_params++] = GPARAM_ASZ;
- rv[num_params++] = GPARAM_BSX;
- rv[num_params++] = GPARAM_BSY;
- rv[num_params++] = GPARAM_BSZ;
- rv[num_params++] = GPARAM_CSX;
- rv[num_params++] = GPARAM_CSY;
- rv[num_params++] = GPARAM_CSZ;
-// rv[num_params++] = GPARAM_R;
-// rv[num_params++] = GPARAM_DIV;
-
- M = gsl_matrix_calloc(num_params, num_params);
- v = gsl_vector_calloc(num_params);
-
- reflections = crystal_get_reflections(cr);
- G = crystal_get_osf(cr);
- B = crystal_get_Bfac(cr);
-
- /* Post-refinement terms */
- for ( refl = first_refl(reflections, &iter);
+
+ n = reflist_new();
+
+ for ( refl = first_refl(input, &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
- signed int ha, ka, la;
- double I_full, delta_I, esd;
- double w;
- double I_partial;
- int k;
- double p, L, s;
- double fx;
- Reflection *match;
- double gradients[num_params];
-
- if ( get_flag(refl) ) continue;
+ signed int h, k, l;
+ Reflection *rn;
- get_indices(refl, &ha, &ka, &la);
- match = find_refl(full, ha, ka, la);
- if ( match == NULL ) continue;
- I_full = get_intensity(match);
+ get_indices(refl, &h, &k, &l);
+ get_equiv(amb, NULL, idx, h, k, l, &h, &k, &l);
+ get_asymm(sym, h, k, l, &h, &k, &l);
+ rn = add_refl(n, h, k, l);
+ copy_data(rn, refl);
+
+ get_symmetric_indices(rn, &h, &k, &l);
+ get_equiv(amb, NULL, idx, h, k, l, &h, &k, &l);
+ set_symmetric_indices(rn, h, k, l);
+ }
- if ( get_redundancy(match) < 2 ) continue;
+ return n;
+}
- p = get_partiality(refl);
- L = get_lorentz(refl);
- I_partial = get_intensity(refl);
- esd = get_esd_intensity(refl);
- s = resolution(crystal_get_cell(cr), ha, ka, la);
- if ( I_partial < 3.0*esd ) continue;
+static void reindex_cell(UnitCell *cell, SymOpList *amb, int idx)
+{
+ signed int h, k, l;
+ struct rvec na, nb, nc;
+ struct rvec as, bs, cs;
+
+ cell_get_reciprocal(cell, &as.u, &as.v, &as.w,
+ &bs.u, &bs.v, &bs.w,
+ &cs.u, &cs.v, &cs.w);
+
+ get_equiv(amb, NULL, idx, 1, 0, 0, &h, &k, &l);
+ na.u = as.u*h + bs.u*k + cs.u*l;
+ na.v = as.v*h + bs.v*k + cs.v*l;
+ na.w = as.w*h + bs.w*k + cs.w*l;
+
+ get_equiv(amb, NULL, idx, 0, 1, 0, &h, &k, &l);
+ nb.u = as.u*h + bs.u*k + cs.u*l;
+ nb.v = as.v*h + bs.v*k + cs.v*l;
+ nb.w = as.w*h + bs.w*k + cs.w*l;
+
+ get_equiv(amb, NULL, idx, 0, 0, 1, &h, &k, &l);
+ nc.u = as.u*h + bs.u*k + cs.u*l;
+ nc.v = as.v*h + bs.v*k + cs.v*l;
+ nc.w = as.w*h + bs.w*k + cs.w*l;
+
+ cell_set_reciprocal(cell, na.u, na.v, na.w,
+ nb.u, nb.v, nb.w,
+ nc.u, nc.v, nc.w);
+}
- /* Calculate the weight for this reflection */
- w = (s/1e9)*(s/1e9) / (esd*esd);
- /* Calculate all gradients for this reflection */
- for ( k=0; k<num_params; k++ ) {
- gradients[k] = gradient(cr, rv[k], refl, pmodel);
- gradients[k] *= exp(G)*exp(-B*s*s)*I_full/L;
- }
+void try_reindex(Crystal *crin, const RefList *full,
+ SymOpList *sym, SymOpList *amb)
+{
+ RefList *list;
+ Crystal *cr;
+ UnitCell *cell;
+ double residual_original, residual_flipped;
+ int idx, n;
- for ( k=0; k<num_params; k++ ) {
+ if ( sym == NULL || amb == NULL ) return;
- int g;
- double v_c, v_curr;
+ residual_original = residual(crin, full, 0, NULL, NULL, 0);
- for ( g=0; g<num_params; g++ ) {
+ cr = crystal_copy(crin);
- double M_c, M_curr;
+ n = num_equivs(amb, NULL);
- /* Matrix is symmetric */
- if ( g > k ) continue;
+ for ( idx=0; idx<n; idx++ ) {
- M_c = w * gradients[g] * gradients[k];
+ cell = cell_new_from_cell(crystal_get_cell(crin));
+ if ( cell == NULL ) return;
+ reindex_cell(cell, amb, idx);
+ crystal_set_cell(cr, cell);
- M_curr = gsl_matrix_get(M, k, g);
- gsl_matrix_set(M, k, g, M_curr + M_c);
- gsl_matrix_set(M, g, k, M_curr + M_c);
+ list = reindex_reflections(crystal_get_reflections(crin),
+ amb, sym, idx);
+ crystal_set_reflections(cr, list);
- }
+ update_predictions(cr);
+ calculate_partialities(cr, PMODEL_XSPHERE);
- fx = exp(G)*p*exp(-B*s*s)*I_full/L;
- delta_I = I_partial - fx;
- v_c = w * delta_I * gradients[k];
- v_curr = gsl_vector_get(v, k);
- gsl_vector_set(v, k, v_curr + v_c);
+ residual_flipped = residual(cr, full, 0, NULL, NULL, 1);
+ if ( residual_flipped < residual_original ) {
+ crystal_set_cell(crin, cell);
+ crystal_set_reflections(crin, list);
+ residual_original = residual_flipped;
+ } else {
+ cell_free(crystal_get_cell(cr));
+ reflist_free(crystal_get_reflections(cr));
}
-
- nref++;
}
- if ( verbose ) {
- //STATUS("The original equation:\n");
- //show_matrix_eqn(M, v);
- STATUS("%i reflections went into the equations.\n", nref);
- }
-
- if ( nref < num_params ) {
- crystal_set_user_flag(cr, PRFLAG_FEWREFL);
- gsl_matrix_free(M);
- gsl_vector_free(v);
- return 0.0;
- }
-
- max_shift = 0.0;
- shifts = solve_svd(v, M, n_filtered, 0);
- if ( shifts != NULL ) {
-
- for ( param=0; param<num_params; param++ ) {
- double shift = gsl_vector_get(shifts, param);
- if ( verbose ) STATUS("Shift %i: %e\n", param, shift);
- if ( isnan(shift) ) {
- //ERROR("NaN shift parameter %i (image ser %i), "
- // "%i reflections.\n", rv[param],
- // crystal_get_image(cr)->serial,
- // nref);
- } else {
- apply_shift(cr, rv[param], shift);
- }
- if ( fabs(shift) > max_shift ) max_shift = fabs(shift);
- }
+ crystal_free(cr);
+}
+
+
+void write_test_logs(Crystal *crystal, const RefList *full,
+ signed int cycle, int serial)
+{
+ FILE *fh;
+ struct image *image = crystal_get_image(crystal);
+ char tmp[256];
+ char ins[5];
+
+ snprintf(tmp, 256, "pr-logs/parameters-crystal%i.dat", serial);
+
+ if ( cycle == 0 ) {
+ fh = fopen(tmp, "w");
} else {
- crystal_set_user_flag(cr, PRFLAG_SOLVEFAIL);
+ fh = fopen(tmp, "a");
+ }
+
+ if ( fh == NULL ) {
+ ERROR("Failed to open '%s'\n", tmp);
+ return;
}
- gsl_matrix_free(M);
- gsl_vector_free(v);
- gsl_vector_free(shifts);
+ if ( cycle == 0 ) {
+ fprintf(fh, "Image: %s %s\n",
+ image->filename, get_event_string(image->event));
+ }
- return max_shift;
+ if ( cycle >= 0 ) {
+ snprintf(ins, 4, "%i", cycle);
+ } else {
+ ins[0] = 'F';
+ ins[1] = '\0';
+ }
+
+ fprintf(fh, "%s rlp_size = %e\n", ins, crystal_get_profile_radius(crystal)/1e10);
+ fprintf(fh, "%s mosaicity = %e\n", ins, crystal_get_mosaicity(crystal));
+ fprintf(fh, "%s wavelength = %f A\n", ins, image->lambda*1e10);
+ fprintf(fh, "%s bandwidth = %f\n", ins, image->bw);
+ fprintf(fh, "%s my scale factor = %e\n", ins, crystal_get_osf(crystal));
+
+ double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz;
+ cell_get_reciprocal(crystal_get_cell(crystal), &asx, &asy, &asz,
+ &bsx, &bsy, &bsz,
+ &csx, &csy, &csz);
+ fprintf(fh, "%s %f, %f, %f, %f, %f, %f, %f, %f, %f\n", ins,
+ asx/1e10, bsx/1e10, csx/1e10,
+ asy/1e10, bsy/1e10, csy/1e10,
+ asz/1e10, bsz/1e10, csz/1e10);
+ fclose(fh);
}
-double residual(Crystal *cr, const RefList *full, int free,
- int *pn_used, const char *filename)
+void write_specgraph(Crystal *crystal, const RefList *full,
+ signed int cycle, int serial)
{
- double dev = 0.0;
- double G, B;
+ FILE *fh;
+ char tmp[256];
Reflection *refl;
RefListIterator *iter;
- FILE *fh = NULL;
- int n_used = 0;
+ double G = crystal_get_osf(crystal);
+ double B = crystal_get_Bfac(crystal);
+ UnitCell *cell;
+ struct image *image = crystal_get_image(crystal);
+ char ins[5];
- if ( filename != NULL ) {
- fh = fopen(filename, "a");
- if ( fh == NULL ) {
- ERROR("Failed to open '%s'\n", filename);
- }
+ snprintf(tmp, 256, "pr-logs/specgraph-crystal%i.dat", serial);
+
+ if ( cycle == 0 ) {
+ fh = fopen(tmp, "w");
+ } else {
+ fh = fopen(tmp, "a");
}
- G = crystal_get_osf(cr);
- B = crystal_get_Bfac(cr);
+ if ( fh == NULL ) {
+ ERROR("Failed to open '%s'\n", tmp);
+ return;
+ }
- for ( refl = first_refl(crystal_get_reflections(cr), &iter);
+ if ( cycle == 0 ) {
+ fprintf(fh, "Image: %s %s\n",
+ image->filename, get_event_string(image->event));
+ fprintf(fh, "khalf/m 1/d(m) pcalc pobs iteration h k l\n");
+ }
+
+ cell = crystal_get_cell(crystal);
+
+ if ( cycle >= 0 ) {
+ snprintf(ins, 4, "%i", cycle);
+ } else {
+ ins[0] = 'F';
+ ins[1] = '\0';
+ }
+
+ for ( refl = first_refl(crystal_get_reflections(crystal), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
- double p, L, s, w;
+ double corr, Ipart, Ifull, pobs, pcalc;
+ double res;
signed int h, k, l;
Reflection *match;
- double esd, I_full, I_partial;
- double fx, dc;
-
- if ( free && !get_flag(refl) ) continue;
get_indices(refl, &h, &k, &l);
+ res = resolution(cell, h, k, l);
+
match = find_refl(full, h, k, l);
if ( match == NULL ) continue;
- I_full = get_intensity(match);
- if ( get_redundancy(match) < 2 ) continue;
+ /* Don't calculate pobs if reference reflection is weak */
+ if ( fabs(get_intensity(match)) / get_esd_intensity(match) < 3.0 ) continue;
- p = get_partiality(refl);
- L = get_lorentz(refl);
- I_partial = get_intensity(refl);
- esd = get_esd_intensity(refl);
- s = resolution(crystal_get_cell(cr), h, k, l);
-
- if ( I_partial < 3.0*esd ) continue;
+ corr = G * exp(B*res*res) * get_lorentz(refl);
+ Ipart = get_intensity(refl) * corr;
+ Ifull = get_intensity(match);
+ pobs = Ipart / Ifull;
+ pcalc = get_partiality(refl);
- fx = exp(G)*p*exp(-B*s*s)*I_full/L;
- dc = I_partial - fx;
- w = (s/1e9)*(s/1e9)/(esd*esd);
- dev += w*dc*dc;
- n_used++;
+ fprintf(fh, "%e %e %f %f %s %4i %4i %4i\n",
+ get_khalf(refl), 2.0*res, pcalc, pobs, ins, h, k, l);
- if ( fh != NULL ) {
- fprintf(fh, "%4i %4i %4i %e %e\n",
- h, k, l, s, dev);
- }
}
- if ( fh != NULL ) fclose(fh);
-
- if ( pn_used != NULL ) *pn_used = n_used;
- return dev;
+ fclose(fh);
}
-static void write_residual_graph(Crystal *cr, const RefList *full)
+static gsl_multimin_fminimizer *setup_minimiser(Crystal *cr, const RefList *full,
+ int verbose, int serial,
+ struct rf_priv *priv)
{
+ gsl_multimin_fminimizer *min;
+ int n_params = 0;
int i;
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
- double a;
- double step;
- double orig_asx;
- FILE *fh;
- UnitCell *cell;
- cell = crystal_get_cell(cr);
+ /* The parameters to be refined */
+ priv->rv[n_params++] = GPARAM_ANG1;
+ priv->rv[n_params++] = GPARAM_ANG2;
+ priv->rv[n_params++] = GPARAM_R;
+ priv->rv[n_params++] = GPARAM_WAVELENGTH;
+ priv->rv[n_params] = GPARAM_EOL; /* End of list */
+
+ priv->cr = cr;
+ priv->full = full;
+ priv->verbose = verbose;
+ priv->serial = serial;
+
+ priv->f.f = residual_f;
+ priv->f.n = n_params;
+ priv->f.params = priv;
+
+ priv->initial = gsl_vector_alloc(n_params);
+ priv->vals = gsl_vector_alloc(n_params);
+ priv->step = gsl_vector_alloc(n_params);
+
+ for ( i=0; i<n_params; i++ ) {
+ gsl_vector_set(priv->initial, i, get_initial_param(cr, priv->rv[i]));
+ gsl_vector_set(priv->vals, i, 0.0);
+ gsl_vector_set(priv->step, i, 1.0);
+ }
- fh = fopen("residual-plot.dat", "w");
+ min = gsl_multimin_fminimizer_alloc(gsl_multimin_fminimizer_nmsimplex2,
+ n_params);
+ gsl_multimin_fminimizer_set(min, &priv->f, priv->vals, priv->step);
- cell_get_reciprocal(cell, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz);
- a = modulus(asx, asy, asz);
- orig_asx = asx;
+ return min;
+}
- step = a/100.0;
- for ( i=-50; i<=50; i++ ) {
+static void write_grid(Crystal *cr, const RefList *full,
+ signed int cycle, int serial,
+ const enum gparam par1, const enum gparam par2,
+ const char *name)
+{
+ FILE *fh;
+ char fn[64];
+ char ins[5];
+ gsl_multimin_fminimizer *min;
+ struct rf_priv priv;
+ int idx1, idx2;
+ int i;
- double res;
- int n;
- asx = orig_asx + (i*step);
- cell_set_reciprocal(cell, asx, asy, asz,
- bsx, bsy, bsz,
- csx, csy, csz);
- update_predictions(cr);
- calculate_partialities(cr, PMODEL_SCSPHERE);
- res = residual(cr, full, 0, &n, NULL);
- fprintf(fh, "%i %e %e %i\n", i, asx, res, n);
+ min = setup_minimiser(cr, full, 0, serial, &priv);
+
+ idx1 = 99;
+ idx2 = 99;
+ for ( i=0; i<priv.f.n; i++ ) {
+ if ( priv.rv[i] == par1 ) idx1 = i;
+ if ( priv.rv[i] == par2 ) idx2 = i;
}
+ assert(idx1 != 99);
+ assert(idx2 != 99);
- cell_set_reciprocal(cell, orig_asx, asy, asz,
- bsx, bsy, bsz,
- csx, csy, csz);
- update_predictions(cr);
- calculate_partialities(cr, PMODEL_SCSPHERE);
+ if ( cycle >= 0 ) {
+ snprintf(ins, 4, "%i", cycle);
+ } else {
+ ins[0] = 'F';
+ ins[1] = '\0';
+ }
+
+ snprintf(fn, 63, "pr-logs/grid-crystal%i-cycle%s-%s.dat",
+ serial, ins, name);
+ fh = fopen(fn, "w");
+ if ( fh != NULL ) {
+ double v1, v2;
+ fprintf(fh, "%e %e %e %s\n",
+ -5.0*get_scale(par1)+get_initial_param(cr, par1),
+ 5.0*get_scale(par1)+get_initial_param(cr, par1),
+ get_initial_param(cr, par1),
+ get_label(par1));
+ fprintf(fh, "%e %e %e %s\n",
+ -5.0*get_scale(par2)+get_initial_param(cr, par2),
+ 5.0*get_scale(par2)+get_initial_param(cr, par2),
+ get_initial_param(cr, par2),
+ get_label(par2));
+ for ( v2=-5.0; v2<=5.0; v2+=0.25 ) {
+ int first=1;
+ for ( v1=-5.0; v1<=5.0; v1+=0.25 ) {
+ double res;
+ gsl_vector_set(min->x, idx1, v1);
+ gsl_vector_set(min->x, idx2, v2);
+ res = residual_f(min->x, &priv);
+ if ( !first ) fprintf(fh, " ");
+ first = 0;
+ fprintf(fh, "%e", res);
+ }
+ fprintf(fh, "\n");
+ }
+ }
fclose(fh);
+
+ gsl_multimin_fminimizer_free(min);
+ gsl_vector_free(priv.initial);
+ gsl_vector_free(priv.vals);
+ gsl_vector_free(priv.step);
+}
+
+
+void write_gridscan(Crystal *cr, const RefList *full,
+ signed int cycle, int serial)
+{
+ write_grid(cr, full, cycle, serial,
+ GPARAM_ANG1, GPARAM_ANG2, "ang1-ang2");
+ write_grid(cr, full, cycle, serial,
+ GPARAM_ANG1, GPARAM_WAVELENGTH, "ang1-wave");
+ write_grid(cr, full, cycle, serial,
+ GPARAM_R, GPARAM_WAVELENGTH, "R-wave");
}
static void do_pr_refine(Crystal *cr, const RefList *full,
- PartialityModel pmodel, int verbose)
+ PartialityModel pmodel, int verbose, int serial,
+ int cycle, int write_logs,
+ SymOpList *sym, SymOpList *amb)
{
- int i, done;
- double old_dev;
- UnitCell *cell = crystal_get_cell(cr);
+ gsl_multimin_fminimizer *min;
+ struct rf_priv priv;
+ int n_iter = 0;
+ int status;
+ int r;
+ double G;
+ double residual_start, residual_free_start;
+ FILE *fh = NULL;
+
+ try_reindex(cr, full, sym, amb);
- old_dev = residual(cr, full, 0, 0, NULL);
+ residual_start = residual(cr, full, 0, NULL, NULL, 0);
+ residual_free_start = residual(cr, full, 1, NULL, NULL, 0);
if ( verbose ) {
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
- cell_get_reciprocal(cell, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz);
- STATUS("Initial asx = %e\n", asx);
- STATUS("PR initial dev = %10.5e, free dev = %10.5e\n",
- old_dev, residual(cr, full, 1, NULL, NULL));
- }
-
- i = 0;
- done = 0;
- do {
+ STATUS("\nPR initial: dev = %10.5e, free dev = %10.5e\n",
+ residual_start, residual_free_start);
+ }
- double dev;
+ min = setup_minimiser(cr, full, verbose, serial, &priv);
- pr_iterate(cr, full, pmodel, NULL, verbose);
+ if ( verbose ) {
+ double res = residual_f(min->x, &priv);
+ double size = gsl_multimin_fminimizer_size(min);
+ STATUS("At start: %f %f %f %f ----> %f %f %e %f residual = %e size %f\n",
+ gsl_vector_get(min->x, 0),
+ gsl_vector_get(min->x, 1),
+ gsl_vector_get(min->x, 2),
+ gsl_vector_get(min->x, 3),
+ rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)),
+ rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)),
+ get_actual_val(min->x, priv.initial, priv.rv, 2),
+ get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10,
+ res, size);
+ }
- update_predictions(cr);
- calculate_partialities(cr, pmodel);
+ if ( write_logs ) {
- dev = residual(cr, full, 0, 0, NULL);
- if ( fabs(dev - old_dev) < dev*0.0001 ) done = 1;
+ char fn[64];
+
+ snprintf(fn, 63, "pr-logs/crystal%i-cycle%i.log", serial, cycle);
+ fh = fopen(fn, "w");
+ if ( fh != NULL ) {
+ fprintf(fh, "iteration RtoReference CCtoReference nref "
+ "ang1 ang2 radius wavelength");
+ double res = residual_f(min->x, &priv);
+ fprintf(fh, "%5i %10.8f %10.8f %5i %10.8f %10.8f %e %e\n",
+ n_iter, res, 0.0, 0,
+ rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)),
+ rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)),
+ get_actual_val(min->x, priv.initial, priv.rv, 2),
+ get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10);
+ }
+
+ }
+
+ do {
+ double res;
+
+ n_iter++;
+
+ status = gsl_multimin_fminimizer_iterate(min);
+ if ( status ) break;
+
+ res = residual_f(min->x, &priv);
+ if ( isnan(res) ) {
+ status = GSL_ENOPROG;
+ break;
+ }
if ( verbose ) {
- STATUS("PR iter %2i: dev = %10.5e, free dev = %10.5e\n",
- i+1, dev, residual(cr, full, 1, NULL, NULL));
+ double res = residual_f(min->x, &priv);
+ double size = gsl_multimin_fminimizer_size(min);
+ STATUS("%f %f %f %f ----> %f %f %e %f residual = %e size %f\n",
+ gsl_vector_get(min->x, 0),
+ gsl_vector_get(min->x, 1),
+ gsl_vector_get(min->x, 2),
+ gsl_vector_get(min->x, 3),
+ rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)),
+ rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)),
+ get_actual_val(min->x, priv.initial, priv.rv, 2),
+ get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10,
+ res, size);
}
- i++;
- old_dev = dev;
+ if ( fh != NULL ) {
+ fprintf(fh, "%5i %10.8f %10.8f %5i %10.8f %10.8f %e %e\n",
+ n_iter, res, 0.0, 0,
+ rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)),
+ rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)),
+ get_actual_val(min->x, priv.initial, priv.rv, 2),
+ get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10);
+ }
- } while ( i < 30 && !done );
+ status = gsl_multimin_test_size(min->size, 0.005);
+
+ } while ( status == GSL_CONTINUE && n_iter < 1000 );
if ( verbose ) {
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
- cell_get_reciprocal(cell, &asx, &asy, &asz,
- &bsx, &bsy, &bsz,
- &csx, &csy, &csz);
- STATUS("Final asx = %e\n", asx);
+ STATUS("Done with refinement after %i iter\n", n_iter);
+ STATUS("status = %i (%s)\n", status, gsl_strerror(status));
}
-}
+ if ( status == GSL_SUCCESS ) {
-static struct prdata pr_refine(Crystal *cr, const RefList *full,
- PartialityModel pmodel)
-{
- int verbose = 0;
- struct prdata prdata;
+ if ( check_angle_shifts(min->x, priv.initial, priv.rv, &priv) ) return;
- prdata.refined = 0;
- prdata.n_filtered = 0;
+ if ( verbose ) {
- if ( verbose ) {
- write_residual_graph(cr, full);
+ double res = residual_f(min->x, &priv);
+ double size = gsl_multimin_fminimizer_size(min);
+ STATUS("At end: %f %f %f %f ----> %f %f %e %f residual = %e size %f\n",
+ gsl_vector_get(min->x, 0),
+ gsl_vector_get(min->x, 1),
+ gsl_vector_get(min->x, 2),
+ gsl_vector_get(min->x, 3),
+ rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)),
+ rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)),
+ get_actual_val(min->x, priv.initial, priv.rv, 2),
+ get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10,
+ res, size);
+
+ }
+
+ if ( fh != NULL ) {
+ double res = residual_f(min->x, &priv);
+ fprintf(fh, "%5i %10.8f %10.8f %5i %10.8f %10.8f %e %e\n",
+ n_iter, res, 0.0, 0,
+ rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)),
+ rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)),
+ get_actual_val(min->x, priv.initial, priv.rv, 2),
+ get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10);
+ fclose(fh);
+ }
+
+ /* Apply the final shifts */
+ apply_parameters(min->x, priv.initial, priv.rv, cr);
+ update_predictions(cr);
+ calculate_partialities(cr, PMODEL_XSPHERE);
+ r = linear_scale(full, crystal_get_reflections(cr), &G, 0);
+ if ( r == 0 ) {
+ crystal_set_osf(cr, G);
+ }
+
+ if ( verbose ) {
+
+ STATUS("After applying final shifts:\n");
+ STATUS("PR final: dev = %10.5e, free dev = %10.5e\n",
+ residual(cr, full, 0, NULL, NULL, 0),
+ residual(cr, full, 1, NULL, NULL, 0));
+ STATUS("Final R = %e m^-1\n", crystal_get_profile_radius(cr));
+
+ }
+
+ } else {
+ ERROR("Bad refinement: crystal %i\n", serial);
}
- do_pr_refine(cr, full, pmodel, verbose);
+ if ( write_logs ) {
+ write_gridscan(cr, full, cycle, serial);
+ write_specgraph(cr, full, cycle, serial);
+ write_test_logs(cr, full, cycle, serial);
+ }
- if ( crystal_get_user_flag(cr) == 0 ) {
- prdata.refined = 1;
+ if ( crystal_get_profile_radius(cr) > 5e9 ) {
+ ERROR("Very large radius: crystal %i\n", serial);
}
- return prdata;
+ gsl_multimin_fminimizer_free(min);
+ gsl_vector_free(priv.initial);
+ gsl_vector_free(priv.vals);
+ gsl_vector_free(priv.step);
}
@@ -706,7 +943,13 @@ struct refine_args
RefList *full;
Crystal *crystal;
PartialityModel pmodel;
+ int serial;
struct prdata prdata;
+ int verbose;
+ int cycle;
+ int no_logs;
+ SymOpList *sym;
+ SymOpList *amb;
};
@@ -724,8 +967,18 @@ static void refine_image(void *task, int id)
{
struct refine_args *pargs = task;
Crystal *cr = pargs->crystal;
+ int write_logs = 0;
+
+ write_logs = !pargs->no_logs && (pargs->serial % 20 == 0);
+ pargs->prdata.refined = 0;
- pargs->prdata = pr_refine(cr, pargs->full, pargs->pmodel);
+ do_pr_refine(cr, pargs->full, pargs->pmodel, pargs->verbose,
+ pargs->serial, pargs->cycle, write_logs,
+ pargs->sym, pargs->amb);
+
+ if ( crystal_get_user_flag(cr) == 0 ) {
+ pargs->prdata.refined = 1;
+ }
}
@@ -738,6 +991,7 @@ static void *get_image(void *vqargs)
memcpy(task, &qargs->task_defaults, sizeof(struct refine_args));
task->crystal = qargs->crystals[qargs->n_started];
+ task->serial = qargs->n_started;
qargs->n_started++;
@@ -757,7 +1011,9 @@ static void done_image(void *vqargs, void *task)
void refine_all(Crystal **crystals, int n_crystals,
- RefList *full, int nthreads, PartialityModel pmodel)
+ RefList *full, int nthreads, PartialityModel pmodel,
+ int verbose, int cycle, int no_logs,
+ SymOpList *sym, SymOpList *amb)
{
struct refine_args task_defaults;
struct queue_args qargs;
@@ -766,7 +1022,11 @@ void refine_all(Crystal **crystals, int n_crystals,
task_defaults.crystal = NULL;
task_defaults.pmodel = pmodel;
task_defaults.prdata.refined = 0;
- task_defaults.prdata.n_filtered = 0;
+ task_defaults.verbose = verbose;
+ task_defaults.cycle = cycle;
+ task_defaults.no_logs = no_logs;
+ task_defaults.sym = sym;
+ task_defaults.amb = amb;
qargs.task_defaults = task_defaults;
qargs.n_started = 0;
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 7f395013..e2416a23 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -3,11 +3,11 @@
*
* Post refinement
*
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2015 Thomas White <taw@physics.org>
+ * 2010-2018 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -41,6 +41,7 @@
#include "utils.h"
#include "crystal.h"
#include "geometry.h"
+#include "symmetry.h"
enum prflag
@@ -51,19 +52,31 @@ enum prflag
PRFLAG_EARLY = 18,
PRFLAG_CC = 19,
PRFLAG_BIGB = 20,
+ PRFLAG_SCALEBAD = 21,
};
extern const char *str_prflag(enum prflag flag);
extern void refine_all(Crystal **crystals, int n_crystals,
- RefList *full, int nthreads, PartialityModel pmodel);
+ RefList *full, int nthreads, PartialityModel pmodel,
+ int verbose, int cycle, int no_logs,
+ SymOpList *sym, SymOpList *amb);
+
+extern void write_gridscan(Crystal *cr, const RefList *full,
+ int cycle, int serial);
+
+extern void write_specgraph(Crystal *crystal, const RefList *full,
+ signed int cycle, int serial);
/* Exported so it can be poked by tests/pr_p_gradient_check */
extern double gradient(Crystal *cr, int k, Reflection *refl,
PartialityModel pmodel);
extern double residual(Crystal *cr, const RefList *full, int free,
- int *pn_used, const char *filename);
+ int *pn_used, const char *filename, int complain);
+
+extern void write_test_logs(Crystal *crystal, const RefList *full,
+ signed int cycle, int serial);
#endif /* POST_REFINEMENT_H */
diff --git a/src/process_image.c b/src/process_image.c
index 87089289..dbd4f427 100644
--- a/src/process_image.c
+++ b/src/process_image.c
@@ -246,7 +246,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
if ( iargs->fix_bandwidth >= 0.0 ) {
image.bw = iargs->fix_bandwidth;
} else {
- image.bw = 0.00000001;
+ image.bw = 0.0013;
}
if ( image_feature_count(image.features) < iargs->min_peaks ) {
@@ -301,7 +301,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs,
/* Integrate! */
time_accounts_set(taccs, TACC_INTEGRATION);
sb_shared->pings[cookie]++;
- integrate_all_5(&image, iargs->int_meth, PMODEL_SCSPHERE,
+ integrate_all_5(&image, iargs->int_meth, PMODEL_XSPHERE,
iargs->push_res,
iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
iargs->int_diag, iargs->int_diag_h,
diff --git a/src/rejection.c b/src/rejection.c
index 10f13f9c..9d41f9b4 100644
--- a/src/rejection.c
+++ b/src/rejection.c
@@ -137,7 +137,7 @@ static void check_cc(Crystal *cr, RefList *full)
pcalc = get_partiality(refl);
/* Observed partiality */
- corr = exp(-G) * exp(B*res*res) * get_lorentz(refl);
+ corr = G * exp(B*res*res) * get_lorentz(refl);
Ipart = get_intensity(refl) * corr;
pobs = Ipart / get_intensity(match);
diff --git a/src/scaling.c b/src/scaling.c
index 77a64384..f4c93c66 100644
--- a/src/scaling.c
+++ b/src/scaling.c
@@ -3,11 +3,11 @@
*
* Scaling
*
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2015 Thomas White <taw@physics.org>
+ * 2010-2017 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -38,16 +38,14 @@
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_blas.h>
+#include <gsl/gsl_fit.h>
#include "merge.h"
#include "post-refinement.h"
#include "symmetry.h"
#include "cell.h"
#include "cell-utils.h"
-
-
-/* Maximum number of iterations of NLSq to do for each image per macrocycle. */
-#define MAX_CYCLES (10)
+#include "scaling.h"
/* Apply the given shift to the 'k'th parameter of 'image'. */
@@ -64,9 +62,9 @@ static void apply_shift(Crystal *cr, int k, double shift)
break;
case GPARAM_OSF :
- t = crystal_get_osf(cr);
+ t = -log(crystal_get_osf(cr));
t += shift;
- crystal_set_osf(cr, t);
+ crystal_set_osf(cr, exp(-t));
break;
default :
@@ -78,8 +76,7 @@ static void apply_shift(Crystal *cr, int k, double shift)
/* Perform one cycle of scaling of 'cr' against 'full' */
-static double scale_iterate(Crystal *cr, const RefList *full,
- PartialityModel pmodel, int *nr)
+static double scale_iterate(Crystal *cr, const RefList *full, int *nr)
{
gsl_matrix *M;
gsl_vector *v;
@@ -182,7 +179,7 @@ static double scale_iterate(Crystal *cr, const RefList *full,
}
- fx = G + log(p) - log(L) - B*s*s + log(I_full);
+ fx = -log(G) + log(p) - log(L) - B*s*s + log(I_full);
delta_I = log(I_partial) - fx;
v_c = w * delta_I * gradients[k];
v_curr = gsl_vector_get(v, k);
@@ -271,7 +268,7 @@ double log_residual(Crystal *cr, const RefList *full, int free,
if ( I_full <= 0 ) continue; /* Because log */
if ( p <= 0.0 ) continue; /* Because of log */
- fx = G + log(p) - log(L) - B*s*s + log(I_full);
+ fx = -log(G) + log(p) - log(L) - B*s*s + log(I_full);
dc = log(I_partial) - fx;
w = 1.0;
dev += w*dc*dc;
@@ -290,8 +287,7 @@ double log_residual(Crystal *cr, const RefList *full, int free,
}
-static void do_scale_refine(Crystal *cr, const RefList *full,
- PartialityModel pmodel, int *nr)
+static void do_scale_refine(Crystal *cr, const RefList *full, int *nr)
{
int i, done;
double old_dev;
@@ -304,7 +300,7 @@ static void do_scale_refine(Crystal *cr, const RefList *full,
double dev;
- scale_iterate(cr, full, pmodel, nr);
+ scale_iterate(cr, full, nr);
dev = log_residual(cr, full, 0, 0, NULL);
if ( fabs(dev - old_dev) < dev*0.01 ) done = 1;
@@ -312,7 +308,7 @@ static void do_scale_refine(Crystal *cr, const RefList *full,
i++;
old_dev = dev;
- } while ( i < MAX_CYCLES && !done );
+ } while ( i < 10 && !done );
}
@@ -320,7 +316,6 @@ struct scale_args
{
RefList *full;
Crystal *crystal;
- PartialityModel pmodel;
int n_reflections;
};
@@ -339,8 +334,7 @@ struct queue_args
static void scale_crystal(void *task, int id)
{
struct scale_args *pargs = task;
- do_scale_refine(pargs->crystal, pargs->full, pargs->pmodel,
- &pargs->n_reflections);
+ do_scale_refine(pargs->crystal, pargs->full, &pargs->n_reflections);
}
@@ -381,14 +375,12 @@ static double total_log_r(Crystal **crystals, int n_crystals, RefList *full,
int n = 0;
for ( i=0; i<n_crystals; i++ ) {
-
double r;
if ( crystal_get_user_flag(crystals[i]) ) continue;
r = log_residual(crystals[i], full, 0, NULL, NULL);
if ( isnan(r) ) continue;
total += r;
n++;
-
}
if ( ninc != NULL ) *ninc = n;
@@ -397,8 +389,7 @@ static double total_log_r(Crystal **crystals, int n_crystals, RefList *full,
/* Perform iterative scaling, all the way to convergence */
-void scale_all(Crystal **crystals, int n_crystals, int nthreads,
- PartialityModel pmodel)
+void scale_all(Crystal **crystals, int n_crystals, int nthreads)
{
struct scale_args task_defaults;
struct queue_args qargs;
@@ -406,7 +397,6 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads,
int niter = 0;
task_defaults.crystal = NULL;
- task_defaults.pmodel = pmodel;
qargs.task_defaults = task_defaults;
qargs.n_crystals = n_crystals;
@@ -422,7 +412,7 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads,
double bef_res;
full = merge_intensities(crystals, n_crystals, nthreads,
- pmodel, 2, INFINITY, 0);
+ 2, INFINITY, 0);
old_res = new_res;
bef_res = total_log_r(crystals, n_crystals, full, NULL);
@@ -456,3 +446,178 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads,
ERROR("Too many iterations - giving up!\n");
}
}
+
+
+static void scale_crystal_linear(void *task, int id)
+{
+ struct scale_args *pargs = task;
+ int r;
+ double G;
+
+ /* Simple iterative algorithm */
+ r = linear_scale(pargs->full, crystal_get_reflections(pargs->crystal), &G, 0);
+ if ( r == 0 ) {
+ crystal_set_osf(pargs->crystal, G);
+ } else {
+ crystal_set_user_flag(pargs->crystal, PRFLAG_SCALEBAD);
+ }
+}
+
+
+/* Calculates G, by which list2 should be multiplied to fit list1 */
+int linear_scale(const RefList *list1, const RefList *list2, double *G,
+ int complain_loudly)
+{
+ const Reflection *refl1;
+ const Reflection *refl2;
+ RefListIterator *iter;
+ int max_n = 256;
+ int n = 0;
+ double *x;
+ double *y;
+ double *w;
+ int r;
+ double cov11;
+ double sumsq;
+ int n_esd1 = 0;
+ int n_esd2 = 0;
+ int n_ih1 = 0;
+ int n_ih2 = 0;
+ int n_nan1 = 0;
+ int n_nan2 = 0;
+ int n_inf1 = 0;
+ int n_inf2 = 0;
+ int n_part = 0;
+ int n_nom = 0;
+
+ x = malloc(max_n*sizeof(double));
+ w = malloc(max_n*sizeof(double));
+ y = malloc(max_n*sizeof(double));
+ if ( (x==NULL) || (y==NULL) || (w==NULL) ) {
+ ERROR("Failed to allocate memory for scaling.\n");
+ return 1;
+ }
+
+ int nb = 0;
+ for ( refl2 = first_refl_const(list2, &iter);
+ refl2 != NULL;
+ refl2 = next_refl_const(refl2, iter) )
+ {
+ signed int h, k, l;
+ double Ih1, Ih2;
+ double esd1, esd2;
+ nb++;
+
+ get_indices(refl2, &h, &k, &l);
+ refl1 = find_refl(list1, h, k, l);
+ if ( refl1 == NULL ) {
+ n_nom++;
+ continue;
+ }
+
+ Ih1 = get_intensity(refl1);
+ Ih2 = get_intensity(refl2);
+ esd1 = get_esd_intensity(refl1);
+ esd2 = get_esd_intensity(refl2);
+
+ /* Problem cases in approximate descending order of severity */
+ if ( isnan(Ih1) ) { n_nan1++; continue; }
+ if ( isinf(Ih1) ) { n_inf1++; continue; }
+ if ( isnan(Ih2) ) { n_nan2++; continue; }
+ if ( isinf(Ih2) ) { n_inf2++; continue; }
+ if ( get_partiality(refl2) < 0.3 ) { n_part++; continue; }
+ //if ( Ih1 <= 0.0 ) { n_ih1++; continue; }
+ if ( Ih2 <= 0.0 ) { n_ih2++; continue; }
+ //if ( Ih1 <= 3.0*esd1 ) { n_esd1++; continue; }
+ if ( Ih2 <= 3.0*esd2 ) { n_esd2++; continue; }
+ //if ( get_redundancy(refl1) < 2 ) continue;
+
+ if ( n == max_n ) {
+ max_n *= 2;
+ x = realloc(x, max_n*sizeof(double));
+ y = realloc(y, max_n*sizeof(double));
+ w = realloc(w, max_n*sizeof(double));
+ if ( (x==NULL) || (y==NULL) || (w==NULL) ) {
+ ERROR("Failed to allocate memory for scaling.\n");
+ return 1;
+ }
+ }
+
+ x[n] = Ih2 / get_partiality(refl2);
+ y[n] = Ih1;
+ w[n] = get_partiality(refl2);
+ n++;
+
+ }
+
+ if ( n < 2 ) {
+ if ( complain_loudly ) {
+ ERROR("Not enough reflections for scaling (had %i, but %i remain)\n", nb, n);
+ if ( n_esd1 ) ERROR("%i reference reflection esd\n", n_esd1);
+ if ( n_esd2 ) ERROR("%i subject reflection esd\n", n_esd2);
+ if ( n_ih1 ) ERROR("%i reference reflection intensity\n", n_ih1);
+ if ( n_ih2 ) ERROR("%i subject reflection intensity\n", n_ih2);
+ if ( n_nan1 ) ERROR("%i reference reflection nan\n", n_nan1);
+ if ( n_nan2 ) ERROR("%i subject reflection nan\n", n_nan2);
+ if ( n_inf1 ) ERROR("%i reference reflection inf\n", n_inf1);
+ if ( n_inf2 ) ERROR("%i subject reflection inf\n", n_inf2);
+ if ( n_part ) ERROR("%i subject reflection partiality\n", n_part);
+ if ( n_nom ) ERROR("%i no match in reference list\n", n_nom);
+ }
+ *G = 1.0;
+ return 1;
+ }
+
+ r = gsl_fit_wmul(x, 1, w, 1, y, 1, n, G, &cov11, &sumsq);
+
+ if ( r ) {
+ ERROR("Scaling failed.\n");
+ return 1;
+ }
+
+ if ( isnan(*G) ) {
+
+ if ( complain_loudly ) {
+ ERROR("Scaling gave NaN (%i pairs)\n", n);
+ if ( n < 10 ) {
+ int i;
+ for ( i=0; i<n; i++ ) {
+ STATUS("%i %e %e %e\n", i, x[i], y[i], w[n]);
+ }
+ }
+ }
+
+ *G = 1.0;
+ return 1;
+ }
+
+ free(x);
+ free(y);
+ free(w);
+
+ return 0;
+}
+
+
+void scale_all_to_reference(Crystal **crystals, int n_crystals,
+ RefList *reference, int nthreads)
+{
+ struct scale_args task_defaults;
+ struct queue_args qargs;
+
+ task_defaults.crystal = NULL;
+
+ qargs.task_defaults = task_defaults;
+ qargs.n_crystals = n_crystals;
+ qargs.crystals = crystals;
+
+ /* Don't have threads which are doing nothing */
+ if ( n_crystals < nthreads ) nthreads = n_crystals;
+
+ qargs.task_defaults.full = reference;
+ qargs.n_started = 0;
+ qargs.n_done = 0;
+ qargs.n_reflections = 0;
+ run_threads(nthreads, scale_crystal_linear, get_crystal, done_crystal,
+ &qargs, n_crystals, 0, 0, 0);
+}
diff --git a/src/scaling.h b/src/scaling.h
index 534043d9..514bee03 100644
--- a/src/scaling.h
+++ b/src/scaling.h
@@ -3,11 +3,11 @@
*
* Scaling
*
- * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2015 Thomas White <taw@physics.org>
+ * 2010-2017 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -41,7 +41,12 @@
extern double log_residual(Crystal *cr, const RefList *full, int free,
int *pn_used, const char *filename);
-extern void scale_all(Crystal **crystals, int n_crystals, int nthreads,
- PartialityModel pmodel);
+extern int linear_scale(const RefList *list1, const RefList *list2, double *G,
+ int complain_loudly);
+
+extern void scale_all(Crystal **crystals, int n_crystals, int nthreads);
+
+extern void scale_all_to_reference(Crystal **crystals, int n_crystals,
+ RefList *reference, int nthreads);
#endif /* SCALING_H */
diff --git a/tests/.gitignore b/tests/.gitignore
index 2afe2ec8..083e0eaa 100644
--- a/tests/.gitignore
+++ b/tests/.gitignore
@@ -13,3 +13,4 @@ ring_check
prof2d_check
ambi_check
prediction_gradient_check
+linear_scale_check
diff --git a/tests/linear_scale_check.c b/tests/linear_scale_check.c
new file mode 100644
index 00000000..4cb6a1ba
--- /dev/null
+++ b/tests/linear_scale_check.c
@@ -0,0 +1,77 @@
+/*
+ * linear_scale_check.c
+ *
+ * Check that linear scaling works
+ *
+ * Copyright © 2017 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
+ *
+ * Authors:
+ * 2017 Thomas White <taw@physics.org>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+
+#include <stdlib.h>
+#include <stdio.h>
+
+#include <reflist.h>
+
+#include "../src/scaling.h"
+
+int main(int argc, char *argv[])
+{
+ int fail = 0;
+ int i;
+ gsl_rng *rng;
+ RefList *list1;
+ RefList *list2;
+ double G;
+ int r;
+
+ list1 = reflist_new();
+ list2 = reflist_new();
+
+ rng = gsl_rng_alloc(gsl_rng_mt19937);
+
+ for ( i=0; i<50; i++ ) {
+ signed int h, k, l;
+ Reflection *refl1;
+ Reflection *refl2;
+ double intens;
+ h = gsl_rng_uniform_int(rng, 20) - gsl_rng_uniform_int(rng, 40);
+ k = gsl_rng_uniform_int(rng, 20) - gsl_rng_uniform_int(rng, 40);
+ l = gsl_rng_uniform_int(rng, 20) - gsl_rng_uniform_int(rng, 40);
+ refl1 = add_refl(list1, h, k, l);
+ refl2 = add_refl(list2, h, k, l);
+ intens = gsl_rng_uniform(rng); /* [0,1) */
+ set_intensity(refl1, intens);
+ set_partiality(refl1, 1.0);
+ set_intensity(refl2, intens*2.0);
+ set_partiality(refl2, 1.0);
+ }
+
+ r = linear_scale(list1, list2, &G, 1);
+ STATUS("Scaling result: %i, G = %f\n", r, G);
+
+ return fail;
+}
diff --git a/tests/pr_p_gradient_check.c b/tests/pr_p_gradient_check.c
deleted file mode 100644
index 5322fcca..00000000
--- a/tests/pr_p_gradient_check.c
+++ /dev/null
@@ -1,529 +0,0 @@
-/*
- * pr_p_gradient_check.c
- *
- * Check partiality gradients for post refinement
- *
- * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
- *
- * Authors:
- * 2012-2014 Thomas White <taw@physics.org>
- *
- * This file is part of CrystFEL.
- *
- * CrystFEL is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * CrystFEL is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>.
- *
- */
-
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-
-
-#include <stdlib.h>
-#include <stdio.h>
-#include <gsl/gsl_statistics.h>
-#include <getopt.h>
-
-#include <image.h>
-#include <cell.h>
-#include <cell-utils.h>
-#include <geometry.h>
-#include <reflist.h>
-#include "../src/post-refinement.h"
-
-
-static void scan_partialities(RefList *reflections, RefList *compare,
- int *valid, long double *vals[3], int idx,
- PartialityModel pmodel)
-{
- int i;
- Reflection *refl;
- RefListIterator *iter;
-
- i = 0;
- for ( refl = first_refl(reflections, &iter);
- refl != NULL;
- refl = next_refl(refl, iter) )
- {
- signed int h, k, l;
- Reflection *refl2;
- double rlow, rhigh, p;
-
- get_indices(refl, &h, &k, &l);
- refl2 = find_refl(compare, h, k, l);
- if ( refl2 == NULL ) {
- valid[i] = 0;
- i++;
- continue;
- }
-
- get_partial(refl2, &rlow, &rhigh, &p);
- vals[idx][i] = p;
- if ( unlikely(p < 0.0) ) {
- ERROR("Negative partiality! %3i %3i %3i %f\n",
- h, k, l, p);
- }
-
- i++;
- }
-}
-
-
-static UnitCell *new_shifted_cell(UnitCell *input, int k, double shift)
-{
- UnitCell *cell;
- double asx, asy, asz;
- double bsx, bsy, bsz;
- double csx, csy, csz;
-
- cell = cell_new();
- cell_get_reciprocal(input, &asx, &asy, &asz, &bsx, &bsy, &bsz,
- &csx, &csy, &csz);
- switch ( k )
- {
- case GPARAM_ASX : asx += shift; break;
- case GPARAM_ASY : asy += shift; break;
- case GPARAM_ASZ : asz += shift; break;
- case GPARAM_BSX : bsx += shift; break;
- case GPARAM_BSY : bsy += shift; break;
- case GPARAM_BSZ : bsz += shift; break;
- case GPARAM_CSX : csx += shift; break;
- case GPARAM_CSY : csy += shift; break;
- case GPARAM_CSZ : csz += shift; break;
- }
- cell_set_reciprocal(cell, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
-
- return cell;
-}
-
-
-static void shift_parameter(struct image *image, int k, double shift)
-{
- switch ( k )
- {
- case GPARAM_DIV : image->div += shift; break;
- }
-}
-
-
-static Crystal *new_shifted_crystal(Crystal *cr, int refine, double incr_val)
-{
- Crystal *cr_new;
- double r;
- UnitCell *cell;
-
- cr_new = crystal_copy(cr);
- if ( cr_new == NULL ) {
- ERROR("Failed to allocate crystal.\n");
- return NULL;
- }
-
- crystal_set_image(cr_new, crystal_get_image(cr));
- r = crystal_get_profile_radius(cr_new);
-
- switch ( refine ) {
-
- case GPARAM_ASX :
- case GPARAM_ASY :
- case GPARAM_ASZ :
- case GPARAM_BSX :
- case GPARAM_BSY :
- case GPARAM_BSZ :
- case GPARAM_CSX :
- case GPARAM_CSY :
- case GPARAM_CSZ :
- cell = new_shifted_cell(crystal_get_cell(cr), refine,
- incr_val);
- crystal_set_cell(cr_new, cell);
- break;
-
- case GPARAM_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;
-
- default :
- ERROR("Can't shift %i\n", refine);
- break;
-
- }
-
- return cr_new;
-}
-
-
-static void calc_either_side(Crystal *cr, double incr_val,
- int *valid, long double *vals[3], int refine,
- PartialityModel pmodel)
-{
- RefList *compare;
- struct image *image = crystal_get_image(cr);
-
- if ( (refine != GPARAM_DIV) ) {
-
- Crystal *cr_new;
-
- /* Crystal properties */
- cr_new = new_shifted_crystal(cr, refine, -incr_val);
- compare = predict_to_res(cr_new, largest_q(image));
- crystal_set_reflections(cr_new, compare);
- calculate_partialities(cr_new, pmodel);
- scan_partialities(crystal_get_reflections(cr), compare, valid,
- vals, 0, pmodel);
- cell_free(crystal_get_cell(cr_new));
- crystal_free(cr_new);
- reflist_free(compare);
-
- cr_new = new_shifted_crystal(cr, refine, +incr_val);
- compare = predict_to_res(cr_new, largest_q(image));
- crystal_set_reflections(cr_new, compare);
- calculate_partialities(cr_new, pmodel);
- scan_partialities(crystal_get_reflections(cr), compare, valid,
- vals, 2, pmodel);
- cell_free(crystal_get_cell(cr_new));
- crystal_free(cr_new);
- reflist_free(compare);
-
- } else {
-
- struct image im_moved;
- Crystal *cr_new = crystal_copy(cr);
- crystal_set_image(cr_new, &im_moved);
-
- /* "Image" properties */
- im_moved = *image;
- shift_parameter(&im_moved, refine, -incr_val);
- compare = predict_to_res(cr_new, largest_q(&im_moved));
- crystal_set_reflections(cr_new, compare);
- calculate_partialities(cr_new, pmodel);
- scan_partialities(crystal_get_reflections(cr), compare,
- valid, vals, 0, pmodel);
- reflist_free(compare);
-
- im_moved = *image;
- shift_parameter(&im_moved, refine, +incr_val);
- compare = predict_to_res(cr_new, largest_q(&im_moved));
- crystal_set_reflections(cr_new, compare);
- calculate_partialities(cr_new, pmodel);
- scan_partialities(crystal_get_reflections(cr), compare,
- valid, vals, 2, pmodel);
- reflist_free(compare);
-
- }
-}
-
-
-static double test_gradients(Crystal *cr, double incr_val, int refine,
- const char *str, const char *file,
- PartialityModel pmodel, int quiet, int plot)
-{
- Reflection *refl;
- RefListIterator *iter;
- long double *vals[3];
- int i;
- int *valid;
- int nref;
- int n_good, n_invalid, n_small, n_nan, n_bad;
- RefList *reflections;
- FILE *fh = NULL;
- int ntot = 0;
- double total = 0.0;
- char tmp[32];
- double *vec1;
- double *vec2;
- int n_line;
- double cc;
-
- reflections = predict_to_res(cr, largest_q(crystal_get_image(cr)));
- crystal_set_reflections(cr, reflections);
-
- nref = num_reflections(reflections);
- if ( nref < 10 ) {
- ERROR("Too few reflections found. Failing test by default.\n");
- return 0.0;
- }
-
- vals[0] = malloc(nref*sizeof(long double));
- vals[1] = malloc(nref*sizeof(long double));
- vals[2] = malloc(nref*sizeof(long double));
- if ( (vals[0] == NULL) || (vals[1] == NULL) || (vals[2] == NULL) ) {
- ERROR("Couldn't allocate memory.\n");
- return 0.0;
- }
-
- valid = malloc(nref*sizeof(int));
- if ( valid == NULL ) {
- ERROR("Couldn't allocate memory.\n");
- return 0.0;
- }
- for ( i=0; i<nref; i++ ) valid[i] = 1;
-
- scan_partialities(reflections, reflections, valid, vals, 1, pmodel);
-
- calc_either_side(cr, incr_val, valid, vals, refine, pmodel);
-
- if ( plot ) {
- snprintf(tmp, 32, "gradient-test-%s.dat", file);
- fh = fopen(tmp, "w");
- }
-
- vec1 = malloc(nref*sizeof(double));
- vec2 = malloc(nref*sizeof(double));
- if ( (vec1 == NULL) || (vec2 == NULL) ) {
- ERROR("Couldn't allocate memory.\n");
- return 0.0;
- }
-
- n_invalid = 0; n_good = 0;
- n_nan = 0; n_small = 0; n_bad = 0; n_line = 0;
- i = 0;
- for ( refl = first_refl(reflections, &iter);
- refl != NULL;
- refl = next_refl(refl, iter) )
- {
-
- long double grad1, grad2, grad;
- double cgrad;
- signed int h, k, l;
-
- get_indices(refl, &h, &k, &l);
-
- if ( !valid[i] ) {
- n_invalid++;
- i++;
- } else {
-
- double r1, r2, p;
-
- 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);
-
- get_partial(refl, &r1, &r2, &p);
-
- if ( isnan(cgrad) ) {
- n_nan++;
- continue;
- }
-
- if ( plot ) {
- fprintf(fh, "%e %Le\n", cgrad, grad);
- }
-
- vec1[n_line] = cgrad;
- vec2[n_line] = grad;
- n_line++;
-
- if ( (fabsl(cgrad) < 5e-8) && (fabsl(grad) < 5e-8) ) {
- n_small++;
- continue;
- }
-
- total += fabsl(cgrad - grad);
- ntot++;
-
- if ( !within_tolerance(grad, cgrad, 5.0)
- || !within_tolerance(cgrad, grad, 5.0) )
- {
-
- 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 {
-
- //STATUS("OK %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_good++;
-
- }
-
- }
-
- }
-
- STATUS("%3s: %3i within 5%%, %3i outside, %3i nan, %3i invalid, "
- "%3i small. ", str, n_good, n_bad, n_nan, n_invalid, n_small);
-
- if ( plot ) {
- fclose(fh);
- }
-
- cc = gsl_stats_correlation(vec1, 1, vec2, 1, n_line);
- STATUS("CC = %+f\n", cc);
- return cc;
-}
-
-
-int main(int argc, char *argv[])
-{
- struct image image;
- const double incr_frac = 1.0/1000000.0;
- double incr_val;
- double ax, ay, az;
- double bx, by, bz;
- double cx, cy, cz;
- UnitCell *cell;
- Crystal *cr;
- struct quaternion orientation;
- int i;
- int fail = 0;
- int quiet = 0;
- int plot = 0;
- int c;
- gsl_rng *rng;
-
- 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;
-
- }
-
- }
-
- image.det = simple_geometry(&image, 1024, 1024);
- image.det->panels[0].res = 13333.3;
- image.det->panels[0].clen = 80e-3;
- image.det->panels[0].coffset = 0.0;
-
- image.lambda = ph_en_to_lambda(eV_to_J(8000.0));
- image.div = 1e-3;
- image.bw = 0.01;
- image.filename = malloc(256);
-
- cr = crystal_new();
- if ( cr == NULL ) {
- ERROR("Failed to allocate crystal.\n");
- return 1;
- }
- crystal_set_mosaicity(cr, 0.0);
- crystal_set_profile_radius(cr, 0.005e9);
- crystal_set_image(cr, &image);
-
- cell = cell_new_from_parameters(10.0e-9, 10.0e-9, 10.0e-9,
- deg2rad(90.0),
- deg2rad(90.0),
- deg2rad(90.0));
-
- rng = gsl_rng_alloc(gsl_rng_mt19937);
-
- for ( i=0; i<2; i++ ) {
-
- UnitCell *rot;
- double val;
- PartialityModel pmodel;
-
- if ( i == 0 ) {
- pmodel = PMODEL_SCSPHERE;
- STATUS("Testing SCSphere model:\n");
- } else if ( i == 1 ) {
- pmodel = PMODEL_SCGAUSSIAN;
- STATUS("Testing SCGaussian model.\n");
- } else {
- ERROR("WTF?\n");
- return 1;
- }
-
- orientation = random_quaternion(rng);
- rot = cell_rotate(cell, orientation);
- crystal_set_cell(cr, rot);
-
- cell_get_reciprocal(rot,
- &ax, &ay, &az, &bx, &by,
- &bz, &cx, &cy, &cz);
-
- incr_val = incr_frac * image.div;
- val = test_gradients(cr, incr_val, GPARAM_DIV, "div", "div",
- pmodel, quiet, plot);
- if ( val < 0.99 ) fail = 1;
-
- incr_val = incr_frac * crystal_get_profile_radius(cr);
- val = test_gradients(cr, incr_val, GPARAM_R, "R", "R", pmodel,
- quiet, plot);
- if ( val < 0.99 ) fail = 1;
-
- incr_val = incr_frac * ax;
- val = test_gradients(cr, incr_val, GPARAM_ASX, "ax*", "x", pmodel,
- quiet, plot);
- if ( val < 0.99 ) fail = 1;
- incr_val = incr_frac * bx;
- val = test_gradients(cr, incr_val, GPARAM_BSX, "bx*", "x", pmodel,
- quiet, plot);
- if ( val < 0.99 ) fail = 1;
- incr_val = incr_frac * cx;
- val = test_gradients(cr, incr_val, GPARAM_CSX, "cx*", "x", pmodel,
- quiet, plot);
- if ( val < 0.99 ) fail = 1;
-
- incr_val = incr_frac * ay;
- val = test_gradients(cr, incr_val, GPARAM_ASY, "ay*", "y", pmodel,
- quiet, plot);
- if ( val < 0.99 ) fail = 1;
- incr_val = incr_frac * by;
- val = test_gradients(cr, incr_val, GPARAM_BSY, "by*", "y", pmodel,
- quiet, plot);
- if ( val < 0.99 ) fail = 1;
- incr_val = incr_frac * cy;
- val = test_gradients(cr, incr_val, GPARAM_CSY, "cy*", "y", pmodel,
- quiet, plot);
- if ( val < 0.99 ) fail = 1;
-
- incr_val = incr_frac * az;
- val = test_gradients(cr, incr_val, GPARAM_ASZ, "az*", "z", pmodel,
- quiet, plot);
- if ( val < 0.99 ) fail = 1;
- incr_val = incr_frac * bz;
- val = test_gradients(cr, incr_val, GPARAM_BSZ, "bz*", "z", pmodel,
- quiet, plot);
- if ( val < 0.99 ) fail = 1;
- incr_val = incr_frac * cz;
- val = test_gradients(cr, incr_val, GPARAM_CSZ, "cz*", "z", pmodel,
- quiet, plot);
- if ( val < 0.99 ) fail = 1;
-
- }
-
- gsl_rng_free(rng);
-
- return fail;
-}
diff --git a/tests/prediction_gradient_check.c b/tests/prediction_gradient_check.c
index 0ffc07ca..85d61ead 100644
--- a/tests/prediction_gradient_check.c
+++ b/tests/prediction_gradient_check.c
@@ -73,7 +73,6 @@ static void scan(RefList *reflections, RefList *compare,
{
signed int h, k, l;
Reflection *refl2;
- double rlow, rhigh, p;
double fs, ss, xh, yh;
struct panel *panel;
@@ -85,7 +84,6 @@ static void scan(RefList *reflections, RefList *compare,
continue;
}
- get_partial(refl2, &rlow, &rhigh, &p);
get_detector_pos(refl2, &fs, &ss);
panel = get_panel(refl2);
twod_mapping(fs, ss, &xh, &yh, panel);
@@ -93,7 +91,7 @@ static void scan(RefList *reflections, RefList *compare,
switch ( checkrxy ) {
case 0 :
- vals[idx][i] = (rlow + rhigh)/2.0;
+ vals[idx][i] = get_exerr(refl2);
break;
case 1 :
@@ -279,8 +277,6 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
i++;
} else {
- double r1, r2, p;
-
grad1 = (vals[1][i] - vals[0][i]) / incr_val;
grad2 = (vals[2][i] - vals[1][i]) / incr_val;
grad = (grad1 + grad2) / 2.0;
@@ -300,18 +296,14 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
if ( checkrxy == 1 ) {
cgrad = x_gradient(refine, refl,
crystal_get_cell(cr),
- &image->det->panels[0],
- crystal_get_image(cr)->lambda);
+ &image->det->panels[0]);
} else {
cgrad = y_gradient(refine, refl,
crystal_get_cell(cr),
- &image->det->panels[0],
- crystal_get_image(cr)->lambda);
+ &image->det->panels[0]);
}
}
- get_partial(refl, &r1, &r2, &p);
-
if ( isnan(cgrad) ) {
n_nan++;
continue;
@@ -338,11 +330,11 @@ static double test_gradients(Crystal *cr, double incr_val, int refine,
{
if ( !quiet ) {
- STATUS("!- %s %3i %3i %3i"
- " %10.2Le %10.2e ratio = %5.2Lf"
- " %10.2e %10.2e\n",
+ STATUS("!- %s %3i %3i %3i "
+ "%10.2Le %10.2e "
+ "ratio = %5.2Lf\n",
str, h, k, l, grad, cgrad,
- cgrad/grad, r1, r2);
+ cgrad/grad);
}
n_bad++;
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-30 13:19:10 +0000