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

3 files changed, 118 insertions, 25 deletions

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;
 }