From 11aa16f82432df051f1e055b28315346d9106188 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Thu, 24 Aug 2023 14:35:48 +0200 Subject: Use built-in Mille writer instead of wrapping C++ version --- libcrystfel/src/crystfel-mille.c | 131 ++++++++++++++++++++++++++++++++++----- libcrystfel/src/crystfel-mille.h | 6 +- libcrystfel/src/predict-refine.c | 6 -- 3 files changed, 118 insertions(+), 25 deletions(-) (limited to 'libcrystfel/src') diff --git a/libcrystfel/src/crystfel-mille.c b/libcrystfel/src/crystfel-mille.c index 036848de..aee391ea 100644 --- a/libcrystfel/src/crystfel-mille.c +++ b/libcrystfel/src/crystfel-mille.c @@ -37,11 +37,6 @@ #include "predict-refine.h" #include "profile.h" -#ifdef HAVE_MILLEPEDE -#include -#endif /* HAVE_MILLEPEDE */ - - int mille_label(int group_serial, enum gparam param) { switch ( param ) { @@ -72,7 +67,84 @@ enum gparam mille_unlabel(int n) } -#ifdef HAVE_MILLEPEDE +struct mille +{ + float *float_arr; + int *int_arr; + int max_entries; + int n; + FILE *fh; +}; + +typedef struct mille Mille; + + +static void mille_add_measurement(Mille *m, + int NLC, float *derLc, + int NGL, float *derGl, int *labels, + float rMeas, float sigma) +{ + int space_required; + int i; + + if ( m == NULL ) return; + + /* Allocate extra space if necessary */ + space_required = m->n + NLC + NGL + 2; + if ( space_required > m->max_entries ) { + + float *new_float_arr; + int *new_int_arr; + int new_max_entries; + + if ( m->max_entries == 0 ) { + new_max_entries = 256; + } else { + new_max_entries = m->max_entries; + } + + while ( new_max_entries < space_required ) { + new_max_entries *= 2; + } + + new_float_arr = realloc(m->float_arr, new_max_entries*sizeof(float)); + new_int_arr = realloc(m->int_arr, new_max_entries*sizeof(int)); + if ( (new_float_arr == NULL) || (new_int_arr == NULL) ) return; + + m->float_arr = new_float_arr; + m->int_arr = new_int_arr; + m->max_entries = new_max_entries; + } + + /* The measurement */ + m->float_arr[m->n] = rMeas; + m->int_arr[m->n] = 0; + m->n++; + + /* Local gradients */ + for ( i=0; ifloat_arr[m->n] = derLc[i]; + m->int_arr[m->n] = i+1; + m->n++; + } + } + + /* The measurement error */ + m->float_arr[m->n] = sigma; + m->int_arr[m->n] = 0; + m->n++; + + /* Global gradients */ + for ( i=0; i 0) ) { + m->float_arr[m->n] = derGl[i]; + m->int_arr[m->n] = labels[i]; + m->n++; + } + } +} + void write_mille(Mille *mille, int n, UnitCell *cell, struct reflpeak *rps, struct image *image, @@ -176,29 +248,58 @@ void write_mille(Mille *mille, int n, UnitCell *cell, } -Mille *crystfel_mille_new(const char *outFileName, - int asBinary, - int writeZero) +Mille *crystfel_mille_new(const char *outFileName) { - return mille_new(outFileName, asBinary, writeZero); + Mille *m; + + m = malloc(sizeof(Mille)); + if ( m == NULL ) return NULL; + + m->max_entries = 0; + m->n = 0; + m->float_arr = NULL; + m->int_arr = NULL; + + m->fh = fopen(outFileName, "wb"); + if ( m->fh == NULL ) { + ERROR("Failed to open Mille file '%s'\n", outFileName); + free(m); + return NULL; + } + + + return m; } void crystfel_mille_free(Mille *m) { - mille_free(m); + if ( m == NULL ) return; + fclose(m->fh); + free(m->float_arr); + free(m->int_arr); + free(m); } void crystfel_mille_delete_last_record(Mille *m) { - mille_delete_last_record(m); + m->n = 0; } void crystfel_mille_write_record(Mille *m) { - mille_write_record(m); -} + float nf = 0.0; + int ni = 0; + int nw = (m->n * 2)+2; + + fwrite(&nw, sizeof(int), 1, m->fh); -#endif /* HAVE_MILLEPEDE */ + fwrite(&nf, sizeof(float), 1, m->fh); + fwrite(m->float_arr, sizeof(float), m->n, m->fh); + + fwrite(&ni, sizeof(int), 1, m->fh); + fwrite(m->int_arr, sizeof(int), m->n, m->fh); + m->n = 0; +} diff --git a/libcrystfel/src/crystfel-mille.h b/libcrystfel/src/crystfel-mille.h index e3f4553e..a4b83815 100644 --- a/libcrystfel/src/crystfel-mille.h +++ b/libcrystfel/src/crystfel-mille.h @@ -31,7 +31,7 @@ #include -typedef void *Mille; +typedef struct mille Mille; #include "cell.h" #include "image.h" @@ -42,9 +42,7 @@ typedef void *Mille; * Detector geometry refinement using Millepede */ -extern Mille *crystfel_mille_new(const char *outFileName, - int asBinary, - int writeZero); +extern Mille *crystfel_mille_new(const char *outFileName); extern void crystfel_mille_free(Mille *m); diff --git a/libcrystfel/src/predict-refine.c b/libcrystfel/src/predict-refine.c index 9ca95bfc..58bf3580 100644 --- a/libcrystfel/src/predict-refine.c +++ b/libcrystfel/src/predict-refine.c @@ -868,13 +868,11 @@ int refine_prediction(struct image *image, Crystal *cr, Mille *mille) pred_residual(rps, n, image->detgeom)); crystal_add_notes(cr, tmp); - #ifdef HAVE_MILLEPEDE if ( mille != NULL ) { profile_start("mille-calc"); write_mille(mille, n, crystal_get_cell(cr), rps, image, Minvs); profile_end("mille-calc"); } - #endif for ( i=0; idetgeom->n_panels; i++ ) { gsl_matrix_free(Minvs[i]); @@ -887,21 +885,17 @@ int refine_prediction(struct image *image, Crystal *cr, Mille *mille) n = pair_peaks(image, cr, NULL, rps); free_rps_noreflist(rps, n); if ( n < 10 ) { - #ifdef HAVE_MILLEPEDE if ( mille != NULL ) { crystfel_mille_delete_last_record(mille); } - #endif return 1; } - #ifdef HAVE_MILLEPEDE if ( mille != NULL ) { profile_start("mille-write"); crystfel_mille_write_record(mille); profile_end("mille-write"); } - #endif return 0; } -- cgit v1.2.3