diff options
-rw-r--r-- | libcrystfel/src/spectrum.c | 211 | ||||
-rw-r--r-- | tests/CMakeLists.txt | 5 | ||||
-rw-r--r-- | tests/spectrum_check.c | 77 |
3 files changed, 280 insertions, 13 deletions
diff --git a/libcrystfel/src/spectrum.c b/libcrystfel/src/spectrum.c index 315202a4..3113c75c 100644 --- a/libcrystfel/src/spectrum.c +++ b/libcrystfel/src/spectrum.c @@ -30,6 +30,9 @@ #include <config.h> #endif +#include <assert.h> +#include <gsl/gsl_sort.h> + #include "spectrum.h" #include "utils.h" @@ -154,11 +157,36 @@ struct gaussian spectrum_get_gaussian(Spectrum *s, int n) */ double spectrum_get_density_at_k(Spectrum *s, double k) { + if ( s->rep == SPEC_HISTOGRAM ) { + int i = 0; + double frac; + if ( k <= s->k[0] ) return 0.0; + if ( k >= s->k[s->n_samples-1] ) return 0.0; + /* k is definitely after the first sample, and definitely + * before the last one */ + while ( (s->k[i] < k) && (i<s->n_samples) ) i++; + assert(i < s->n_samples); + frac = (k - s->k[i-1]) / (s->k[i] - s->k[i-1]); + return s->pdf[i-1] + frac * (s->pdf[i] - s->pdf[i-1]); + } + + if ( s->rep == SPEC_GAUSSIANS ) { + double total = 0.0; + int i; + for ( i=0; i<s->n_gaussians; i++ ) { + double a = s->gaussians[i].height; + double b = s->gaussians[i].kcen; + double c = s->gaussians[i].sigma; + total += a*exp(-(k-b)*(k-b)/(2.0*c*c)); + } + return total; + } + return 0.0; } -static double smallest(double *vals, int n_vals) +static double smallest_in_list(double *vals, int n_vals) { int i; double v = +INFINITY; @@ -169,7 +197,7 @@ static double smallest(double *vals, int n_vals) } -static double largest(double *vals, int n_vals) +static double largest_in_list(double *vals, int n_vals) { int i; double v = -INFINITY; @@ -186,7 +214,7 @@ static double gauss_low(struct gaussian *gauss, int n_gauss) double v = +INFINITY; for ( i=0; i<n_gauss; i++ ) { double gv = gauss[i].kcen - 5.0*gauss[i].sigma; - if ( gv > v ) v = gv; + if ( gv < v ) v = gv; } return v; } @@ -216,8 +244,8 @@ static double gauss_high(struct gaussian *gauss, int n_gauss) void spectrum_get_range(Spectrum *s, double *kmin, double *kmax) { if ( s->rep == SPEC_HISTOGRAM ) { - *kmin = smallest(s->k, s->n_samples); - *kmax = largest(s->k, s->n_samples); + *kmin = smallest_in_list(s->k, s->n_samples); + *kmax = largest_in_list(s->k, s->n_samples); } else { assert(s->rep == SPEC_GAUSSIANS); *kmin = gauss_low(s->gaussians, s->n_gaussians); @@ -226,30 +254,152 @@ void spectrum_get_range(Spectrum *s, double *kmin, double *kmax) } +static signed int cmp_gauss(const void *va, const void *vb) +{ + const struct gaussian *a = va; + const struct gaussian *b = vb; + /* Integral of Gaussian = height * std deviation */ + if ( a->height * a->sigma > b->height * b->sigma ) return +1; + return -1; +} + + +static void normalise_gaussians(struct gaussian *gauss, int n_gauss) +{ + int i; + double total_area = 0.0; + for ( i=0; i<n_gauss; i++ ) { + total_area += gauss[i].height * gauss[i].sigma; + } + total_area *= sqrt(2.0 * M_PI); + for ( i=0; i<n_gauss; i++ ) { + gauss[i].height /= total_area; + } +} + + /** * \param s A \ref Spectrum * \param gs Pointer to array of \ref gaussian structures * \param n_gauss Number of Gaussians in \p gs * * Sets the spectrum in terms of a sum of Gaussians. + * + * The spectral density function will be normalised, so only the relative areas + * under the curves are relevant. + * + * The input array will be copied, so you can safely free it after calling this + * function. */ void spectrum_set_gaussians(Spectrum *s, struct gaussian *gs, int n_gauss) { - /* FIXME: sort them */ + /* Free old contents (if any - may be NULL) */ + free(s->gaussians); + free(s->k); + free(s->pdf); + + s->gaussians = malloc(n_gauss * sizeof(struct gaussian)); + if ( s->gaussians == NULL ) return; + + memcpy(s->gaussians, gs, n_gauss*sizeof(struct gaussian)); + s->n_gaussians = n_gauss; + s->rep = SPEC_GAUSSIANS; + + qsort(s->gaussians, s->n_gaussians, sizeof(struct gaussian), cmp_gauss); + normalise_gaussians(s->gaussians, s->n_gaussians); +} + + +/* Samples must already have been sorted */ +static void normalise_histogram(double *k, double *pdf, int n) +{ + int i; + double total_area = 0.0; + double old_k = k[0]; + double old_pdf = pdf[0]; + for ( i=1; i<n; i++ ) { + total_area += (pdf[i]+old_pdf)*(k[i]-old_k)/2.0; + old_k = k[i]; + old_pdf = pdf[i]; + } + + printf("total area under histogram = %f\n", total_area); + + for ( i=0; i<n; i++ ) { + pdf[i] /= total_area; + } } /** * \param s A \ref Spectrum - * \param kcens Pointer to array of k values in the centres of the bins - * \param heights Pointer to array of spectral density values - * \param nbins Number of bins + * \param kvals Pointer to array of k values + * \param heights Pointer to array of spectral density samples + * \param n Number of samples * - * Sets the spectrum in terms of a histogram. + * Sets the spectrum in terms of samples of the probability density function. + * The spectral density function will be normalised. The spectrum can have a + * non-zero value only for k values in the range [kmin,kmax] (exclusive + * interval), where kmin and kmax are the smallest and largest values in + * \p kvals. + * + * The input arrays will be copied, so you can safely free them after calling + * this function. */ -void spectrum_set_histogram(Spectrum *s, double *kcens, double *heights, - int nbins) +void spectrum_set_histogram(Spectrum *s, double *kvals, double *heights, + int n) { + /* Free old contents (if any - may be NULL) */ + free(s->gaussians); + free(s->k); + free(s->pdf); + + s->k = malloc(n * sizeof(double)); + if ( s->k == NULL ) return; + + s->pdf = malloc(n * sizeof(double)); + if ( s->pdf == NULL ) return; + + memcpy(s->k, kvals, n*sizeof(double)); + memcpy(s->pdf, heights, n*sizeof(double)); + s->n_samples = n; + s->rep = SPEC_HISTOGRAM; + + gsl_sort2(s->k, 1, s->pdf, 1, n); + normalise_histogram(s->k, s->pdf, s->n_samples); +} + + +static int read_esrf_spectrum(FILE *fh, Spectrum *s) +{ + double *k = NULL; + double *samp = NULL; + int n_bins = 0; + int max_bins = 0; + + while ( !feof(fh) ) { + + float energy, weight; + if ( fscanf(fh, "%e %e\n", &energy, &weight) != 2 ) return 1; + + if ( n_bins == max_bins ) { + max_bins += 64; + k = realloc(k, max_bins*sizeof(double)); + samp = realloc(samp, max_bins*sizeof(double)); + if ( (k==NULL) || (samp==NULL) ) return 1; + } + + k[n_bins] = ph_eV_to_k(energy*1000.0); + samp[n_bins] = weight; + n_bins++; + + } + + spectrum_set_histogram(s, k, samp, n_bins); + free(k); + free(samp); + + return 0; } @@ -262,5 +412,40 @@ void spectrum_set_histogram(Spectrum *s, double *kcens, double *heights, */ Spectrum *spectrum_load(const char *filename) { - return NULL; + FILE *fh; + Spectrum *s; + char line[1024]; + + fh = fopen(filename, "r"); + if ( fh == NULL ) return NULL; + + s = spectrum_new(); + if ( s == NULL ) { + fclose(fh); + return NULL; + } + + if ( fgets(line, 1024, fh) != line ) { + ERROR("Failed to read '%s'\n", filename); + spectrum_free(s); + return NULL; + } + + chomp(line); + if ( strcmp(line, "# energy/keV current/A") == 0 ) { + if ( read_esrf_spectrum(fh, s) ) { + ERROR("Failed to read ESRF spectrum from %s\n", + filename); + spectrum_free(s); + return NULL; + } + } else { + ERROR("Spectrum format not recognised: %s\n", filename); + fclose(fh); + spectrum_free(s); + return NULL; + } + + fclose(fh); + return s; } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index af22977e..192913ce 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -87,3 +87,8 @@ add_executable(rational_check rational_check.c) target_include_directories(rational_check PRIVATE ${COMMON_INCLUDES}) target_link_libraries(rational_check ${COMMON_LIBRARIES}) add_test(rational_check rational_check) + +add_executable(spectrum_check spectrum_check.c) +target_include_directories(spectrum_check PRIVATE ${COMMON_INCLUDES}) +target_link_libraries(spectrum_check ${COMMON_LIBRARIES}) +add_test(spectrum_check spectrum_check) diff --git a/tests/spectrum_check.c b/tests/spectrum_check.c new file mode 100644 index 00000000..7d37689d --- /dev/null +++ b/tests/spectrum_check.c @@ -0,0 +1,77 @@ +/* + * spectrum_check.c + * + * Check that Spectrum object works + * + * Copyright © 2019 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2019 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 <stdarg.h> + +#include <spectrum.h> +#include <utils.h> + +static int check_integral(Spectrum *s, int nsamp) +{ + double min, max, step; + int i; + double area = 0.0; + + spectrum_get_range(s, &min, &max); + fprintf(stderr, "bounds: %e %e\n", min, max); + step = (max-min)/nsamp; + for ( i=0; i<=nsamp; i++ ) { + double x = min+i*step; + double y = spectrum_get_density_at_k(s, x); + area += y * step; + } + fprintf(stderr, "Total area + %f\n", area); + if ( area > 1.1 ) return 1; + if ( area < 0.9 ) return 1; + return 0; +} + + +int main(int argc, char *argv[]) +{ + Spectrum *s; + struct gaussian gauss; + int r = 0; + + s = spectrum_new(); + gauss.kcen = ph_eV_to_k(9000); + gauss.sigma = ph_eV_to_k(100); + gauss.height = 1.0; + spectrum_set_gaussians(s, &gauss, 1); + r += check_integral(s, 100); + spectrum_free(s); + + return r; +} |