diff options
author | Thomas White <taw@physics.org> | 2023-08-24 14:35:48 +0200 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2023-08-24 14:35:48 +0200 |
commit | 11aa16f82432df051f1e055b28315346d9106188 (patch) | |
tree | 461f5b9f72f3ea69d07b60c8c69ab61f6d73b15c /libcrystfel | |
parent | 9fca86fd0fba293aadaaf88d5d4fc5660b14ef20 (diff) |
Use built-in Mille writer instead of wrapping C++ version
Diffstat (limited to 'libcrystfel')
-rw-r--r-- | libcrystfel/libcrystfel-config.h.meson.in | 1 | ||||
-rw-r--r-- | libcrystfel/meson.build | 3 | ||||
-rw-r--r-- | libcrystfel/src/crystfel-mille.c | 131 | ||||
-rw-r--r-- | libcrystfel/src/crystfel-mille.h | 6 | ||||
-rw-r--r-- | libcrystfel/src/predict-refine.c | 6 |
5 files changed, 118 insertions, 29 deletions
diff --git a/libcrystfel/libcrystfel-config.h.meson.in b/libcrystfel/libcrystfel-config.h.meson.in index c788c017..302fae0b 100644 --- a/libcrystfel/libcrystfel-config.h.meson.in +++ b/libcrystfel/libcrystfel-config.h.meson.in @@ -11,7 +11,6 @@ #mesondefine HAVE_CLOCK_GETTIME #mesondefine HAVE_HDF5 #mesondefine HAVE_SEEDEE -#mesondefine HAVE_MILLEPEDE #mesondefine HAVE_FORKPTY_PTY_H #mesondefine HAVE_FORKPTY_UTIL_H diff --git a/libcrystfel/meson.build b/libcrystfel/meson.build index 81156b61..23b267c1 100644 --- a/libcrystfel/meson.build +++ b/libcrystfel/meson.build @@ -41,9 +41,6 @@ endif millepededep = dependency('millepede', required: false, fallback: ['millepede', 'millepede_dep']) -if millepededep.found() - conf_data.set10('HAVE_MILLEPEDE', true) -endif xgandalfdep = dependency('xgandalf', required: false, 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 <mille_c_wrap.h> -#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; i<NLC; i++ ) { + if ( derLc[i] != 0.0 ) { + m->float_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<NGL; i++ ) { + if ( (derGl[i] != 0.0) && (labels[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 <gsl/gsl_matrix.h> -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; i<image->detgeom->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; } |