Use built-in Mille writer instead of wrapping C++ version

10 files changed, 121 insertions, 118 deletions

diff --git a/config.h.in b/config.h.in
index f0e8c962..c4ac7ebe 100644
--- a/config.h.in
+++ b/config.h.in
@@ -7,4 +7,3 @@
 #mesondefine HAVE_HDF5
 #mesondefine HAVE_ASAPO
 #mesondefine HAVE_SCHED_SETAFFINITY
-#mesondefine HAVE_MILLEPEDE
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;
 }
diff --git a/src/im-sandbox.c b/src/im-sandbox.c
index 6b68326c..27e33b8f 100644
--- a/src/im-sandbox.c
+++ b/src/im-sandbox.c
@@ -366,13 +366,11 @@ static int run_work(const struct index_args *iargs, Stream *st,
 	}
 
 	mille = NULL;
-	#ifdef HAVE_MILLEPEDE
 	if ( iargs->mille ) {
 		char tmp[64];
 		snprintf(tmp, 63, "mille-data-%i.bin", cookie);
-		mille = crystfel_mille_new(tmp, 1, 0);
+		mille = crystfel_mille_new(tmp);
 	}
-	#endif
 
 	while ( !allDone ) {
 
@@ -547,9 +545,7 @@ static int run_work(const struct index_args *iargs, Stream *st,
 		free(pargs.event);
 	}
 
-	#ifdef HAVE_MILLEPEDE
 	crystfel_mille_free(mille);
-	#endif
 
 	/* These are both no-ops if argument is NULL */
 	im_zmq_shutdown(zmqstuff);
diff --git a/subprojects/packagefiles/millepede/meson.build b/subprojects/packagefiles/millepede/meson.build
index 603c709e..2857205f 100644
--- a/subprojects/packagefiles/millepede/meson.build
+++ b/subprojects/packagefiles/millepede/meson.build
@@ -1,5 +1,5 @@
 # Meson file for Millepede-II
-project('millepede', ['c', 'cpp','fortran'],
+project('millepede', ['c', 'fortran'],
         version: '2.0.0',
         license: 'LGPL2+',
         default_options: ['buildtype=debugoptimized'])
@@ -11,11 +11,7 @@ if zlibdep.found()
   pede_c_args += '-DUSE_ZLIB'
 endif
 
-mille = library('mille', ['Mille.cc', 'mille_c_wrap.cc'],
-                install: true)
-
-millepede_dep = declare_dependency(include_directories: '.',
-                                   link_with: mille)
+millepede_dep = declare_dependency()
 
 executable('pede',
            ['mpdef.f90',
diff --git a/subprojects/packagefiles/millepede/mille_c_wrap.cc b/subprojects/packagefiles/millepede/mille_c_wrap.cc
deleted file mode 100644
index 103754a1..00000000
--- a/subprojects/packagefiles/millepede/mille_c_wrap.cc
+++ /dev/null
@@ -1,47 +0,0 @@
-#include "Mille.h"
-#include "mille_c_wrap.h"
-
-
-extern "C" Mille *mille_new(const char *outFileName,
-                            int asBinary,
-                            int writeZero)
-{
-	return new Mille(outFileName, asBinary, writeZero);
-}
-
-
-extern "C" void mille_free(Mille *m)
-{
-	delete m;
-}
-
-
-extern "C" void mille_add_measurement(Mille *m,
-                                      int NLC, const float *derLc,
-                                      int NGL, const  float *derGl,
-                                      const int *label, float rMeas,
-                                      float sigma)
-{
-	m->mille(NLC, derLc, NGL, derGl, label, rMeas, sigma);
-}
-
-
-extern "C" void mille_add_special(Mille *m,
-                                  int nSpecial,
-                                  const float *floatings,
-                                  const int *integers)
-{
-	m->special(nSpecial, floatings, integers);
-}
-
-
-extern "C" void mille_delete_last_record(Mille *m)
-{
-	m->kill();
-}
-
-
-extern "C" void mille_write_record(Mille *m)
-{
-	m->end();
-}
diff --git a/subprojects/packagefiles/millepede/mille_c_wrap.h b/subprojects/packagefiles/millepede/mille_c_wrap.h
deleted file mode 100644
index 1f0fc841..00000000
--- a/subprojects/packagefiles/millepede/mille_c_wrap.h
+++ /dev/null
@@ -1,30 +0,0 @@
-#ifdef __cplusplus
-extern "C" {
-#else
-typedef void *Mille;
-#endif
-
-
-extern Mille *mille_new(const char *outFileName,
-                        int asBinary,
-                        int writeZero);
-
-extern void mille_add_measurement(Mille *m,
-                                  int NLC, const float *derLc,
-                                  int NGL, const  float *derGl,
-                                  const int *label, float rMeas, float sigma);
-
-extern void mille_add_special(Mille *m,
-                              int nSpecial,
-                              const float *floatings,
-                              const int *integers);
-
-extern void mille_delete_last_record(Mille *m);
-
-extern void mille_write_record(Mille *m);
-
-extern void mille_free(Mille *m);
-
-#ifdef __cplusplus
-} /* extern "C" */
-#endif