Refactor MOSFLM autoindexing to use a PTY

1 files changed, 403 insertions, 125 deletions

diff --git a/src/mosflm.c b/src/mosflm.c
index 69d0ea46..080e93a9 100644
--- a/src/mosflm.c
+++ b/src/mosflm.c
@@ -10,6 +10,30 @@
  *
  */
 
+/* TODO:
+ *
+ * Properly read the newmat file (don't use fscanf-- spaces between numers
+ *  are not guaranteed)
+ *
+ * "Success" is indicated by existence of NEWMAT file written by mosflm.
+ *  Better to interact with mosflm directly in order to somehow verify success.
+ *
+ * Investigate how these keywords affect mosflms behavior:
+ *
+ *  MOSAICITY
+ *  DISPERSION
+ *  DIVERGENCE
+ *  POLARISATION
+ *  POSTREF BEAM
+ *  POSTREF USEBEAM OFF
+ *  PREREFINE ON
+ *  EXTRA ON
+ *  POSTREF ON
+ *
+ *  These did not seem to affect the results by my (Rick's) experience, probably
+ *  because they are only used conjunction with image intensity data, but it's
+ *  worth another look at the documentation.
+ */
 
 #ifdef HAVE_CONFIG_H
 #include <config.h>
@@ -23,6 +47,16 @@
 #include <sys/types.h>
 #include <sys/wait.h>
 #include <unistd.h>
+#include <assert.h>
+#include <fcntl.h>
+#include <errno.h>
+
+#if HAVE_FORKPTY_LINUX
+#include <pty.h>
+#elif HAVE_FORKPTY_BSD
+#include <util.h>
+#endif
+
 
 #include "image.h"
 #include "mosflm.h"
@@ -31,40 +65,53 @@
 #include "peaks.h"
 
 
-/*
+#define MOSFLM_VERBOSE 0
 
-todo
 
-** properly read the newmat file (don't use fscanf-- spaces between numers
-   are not guaranteed)
+typedef enum {
+	MOSFLM_INPUT_NONE,
+	MOSFLM_INPUT_LINE,
+	MOSFLM_INPUT_PROMPT
+} MOSFLMInputType;
 
-** "success" is indicated by existence of NEWMAT file written by mosflm.
-   Better to interact with mosflm directly in order to somehow verify success.
 
-** investigate how these keywords affect mosflms behavior:
+struct mosflm_data {
 
-   MOSAICITY
-   DISPERSION
-   DIVERGENCE
-   POLARISATION
-   POSTREF BEAM
-   POSTREF USEBEAM OFF
-   PREREFINE ON
-   EXTRA ON
-   POSTREF ON
+	/* MOSFLM auto-indexing low-level stuff */
+	int                     pty;
+	pid_t                   pid;
+	char                    *rbuffer;
+	int                     rbufpos;
+	int                     rbuflen;
 
-   These did not seem to affect the results by my (Rick's) experience, probably
-   because they are only used conjunction with image intensity data, but it's
-   worth another look at the documentation.
+	/* MOSFLM high-level stuff */
+	char                    newmatfile[128];
+	char                    imagefile[128];
+	char                    sptfile[128];
+	int                     step;
+	UnitCell                *target_cell;
 
-*/
+};
 
 
+static void mosflm_parseline(const char *line, struct image *image,
+                             struct mosflm_data *dirax)
+{
+	#if MOSFLM_VERBOSE
+	char *copy;
 
-#define MOSFLM_VERBOSE 0
+	copy = strdup(line);
+	for ( i=0; i<strlen(copy); i++ ) {
+		if ( copy[i] == '\r' ) copy[i]='r';
+		if ( copy[i] == '\n' ) copy[i]='\0';
+	}
+	STATUS("MOSFLM: %s\n", copy);
+	free(copy);
+	#endif
+}
 
 
-static int read_newmat(const char * filename, struct image *image)
+static int read_newmat(const char *filename, struct image *image)
 {
 	FILE * fh;
 	float asx, asy, asz;
@@ -73,21 +120,21 @@ static int read_newmat(const char * filename, struct image *image)
 	int n;
 	double c;
 
-	fh = fopen(filename,"r");
-	if (fh == NULL){
-		STATUS("found newmat.\n");
+	fh = fopen(filename, "r");
+	if ( fh == NULL ) {
+		STATUS("Found NEWMAT file.\n");
 		return 1;
 	}
-	n  = fscanf(fh,"%f %f %f\n",&asx,&bsx,&csx);
-	n += fscanf(fh,"%f %f %f\n",&asy,&bsy,&csy);
-	n += fscanf(fh,"%f %f %f\n",&asz,&bsz,&csz);
-	if (n != 9) {
+	n  = fscanf(fh, "%f %f %f\n", &asx, &bsx, &csx);
+	n += fscanf(fh, "%f %f %f\n", &asy, &bsy, &csy);
+	n += fscanf(fh, "%f %f %f\n", &asz, &bsz, &csz);
+	if ( n != 9 ) {
 		STATUS("<9 parameters.\n");
 		return 1;
 	}
 	fclose(fh);
 
-	/* mosflm A matrix is multiplied by lambda, so fix this */
+	/* MOSFLM "A" matrix is multiplied by lambda, so fix this */
 	c = 1/image->lambda;
 
 	image->candidate_cells[0] = cell_new();
@@ -103,8 +150,7 @@ static int read_newmat(const char * filename, struct image *image)
 }
 
 
-/* write .spt file for mosflm */
-/* need to sort mosflm peaks by intensity... */
+/* Need to sort mosflm peaks by intensity... */
 struct sptline {
 	double x; /* x coordinate of peak */
 	double y; /* y coordinate of peak */
@@ -124,23 +170,22 @@ static int compare_vals(const void *ap, const void *bp)
 }
 
 
-static void write_spt(struct image *image)
+/* Write .spt file for mosflm */
+static void write_spt(struct image *image, const char *filename)
 {
 	FILE *fh;
 	int i;
-	char filename[1024];
 	double fclen=67.8;  /* fake camera length in mm */
 	double fpix=0.075;  /* fake pixel size in mm */
 	double pix;
 	double height=100;
 	double sigma=1;
 	int nPeaks = image_feature_count(image->features);
-
-	snprintf(filename, 1023, "xfel-%i.spt", image->id);
+	struct sptline *sptlines;
 
 	fh = fopen(filename, "w");
 	if ( !fh ) {
-		ERROR("Couldn't open temporary file xfel.spt\n");
+		ERROR("Couldn't open temporary file '%s'\n", filename);
 		return;
 	}
 
@@ -148,7 +193,6 @@ static void write_spt(struct image *image)
 	fprintf(fh, "%10d %10d\n", 1, 1);
 	fprintf(fh, "%10.5f %10.5f\n", 0.0, 0.0);
 
-	struct sptline *sptlines;
 	sptlines = malloc(sizeof(struct sptline)*nPeaks);
 
 	for ( i=0; i<nPeaks; i++ ) {
@@ -159,7 +203,7 @@ static void write_spt(struct image *image)
 		if ( f == NULL ) continue;
 
 		struct panel *pan;
-		pan = find_panel(image->det,f->x,f->y);
+		pan = find_panel(image->det, f->x, f->y);
 		if ( pan == NULL ) continue;
 
 		pix = 1000/pan->res; /* pixel size in mm */
@@ -189,119 +233,353 @@ static void write_spt(struct image *image)
 }
 
 
-/* Write a dummy 1x1 pixel image file for mosflm.  Without post refinement,
- * mosflm will ignore this, but it must be present. */
-static void write_img(struct image *image)
+/* Write a dummy 1x1 pixel image file for MOSFLM.  Without post refinement,
+ * MOSFLM will ignore this, but it must be present. */
+static void write_img(struct image *image, const char *filename)
 {
 	FILE *fh;
-	char filename[1024];
-	unsigned short int * intimage;
+	unsigned short int *intimage;
 
 	intimage = malloc(sizeof(unsigned short int));
 	intimage[0] = 1;
 
-	snprintf(filename, 1023, "xfel-%i_001.img", image->id);
-
 	fh = fopen(filename, "w");
 	if ( !fh ) {
 		ERROR("Couldn't open temporary file '%s'\n", filename);
 		return;
 	}
 
-	fprintf(fh,"{\nHEADER_BYTES=512;\n");
-	fprintf(fh,"BYTE_ORDER=little_endian;\n");
-	fprintf(fh,"TYPE=unsigned_short;\n");
-	fprintf(fh,"DIM=2;\n");
-	fprintf(fh,"SIZE1=1;\n");
-	fprintf(fh,"SIZE2=1;\n");
-	fprintf(fh,"}\n");
-	while ( ftell(fh) < 512 ) { fprintf(fh," "); };
-	fwrite(intimage,sizeof(unsigned short int),1,fh);
+	/* Write header */
+	fprintf(fh, "{\nHEADER_BYTES=512;\n");
+	fprintf(fh, "BYTE_ORDER=little_endian;\n");
+	fprintf(fh, "TYPE=unsigned_short;\n");
+	fprintf(fh, "DIM=2;\n");
+	fprintf(fh, "SIZE1=1;\n");
+	fprintf(fh, "SIZE2=1;\n");
+	fprintf(fh, "}\n");
+
+	/* Header padding */
+	while ( ftell(fh) < 512 ) fprintf(fh," ");
+
+	fwrite(fh, sizeof(unsigned short int), 1, fh);
+
 	fclose(fh);
 }
 
 
-void run_mosflm(struct image *image, UnitCell *cell)
+static void mosflm_sendline(const char *line, struct mosflm_data *mosflm)
+{
+	#if MOSFLM_VERBOSE
+	char *copy;
+	int i;
+
+	copy = strdup(line);
+	for ( i=0; i<strlen(copy); i++ ) {
+		if ( copy[i] == '\r' ) copy[i]='\0';
+		if ( copy[i] == '\n' ) copy[i]='\0';
+	}
+	STATUS("To MOSFLM: '%s'\n", copy);
+	free(copy);
+	#endif
+
+	write(mosflm->pty, line, strlen(line));
+}
+
+
+static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm)
 {
-	int i,j;
-	char mos_cmd[1024];
-	char symm[64];
+	char tmp[256];
+	char symm[32];
 	const char *sg;
-	double a,b,c,alpha,beta,gamma;
-	double wavelength; /* angstrom */
-	char newmatfile[128];
-	int fail;
-	pid_t pid;
-	int r;
-
-	write_spt(image);
-	write_img(image); /* dummy image */
-
-	wavelength = image->lambda*1e10;
-	cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma);
-	sg = cell_get_spacegroup(cell);
-	sprintf(newmatfile,"xfel-%i.newmat",image->id);
-
-	/* need to remove white space from spacegroup... */
-	j = 0;
-	for(i = 0; i < strlen(sg);i++)
-	{
-		if (sg[i] != ' ') {
-			symm[j] = sg[i];
-			j++;
+	double wavelength;
+	double a, b, c, alpha, beta, gamma;
+	int i, j;
+
+	switch ( mosflm->step ) {
+
+	case 1 :
+		mosflm_sendline("DETECTOR ROTATION HORIZONTAL"
+		                " ANTICLOCKWISE ORIGIN LL FAST HORIZONTAL"
+		                " RECTANGULAR\n", mosflm);
+		break;
+
+	case 2 :
+		cell_get_parameters(mosflm->target_cell, &a, &b, &c,
+		                    &alpha, &beta, &gamma);
+		snprintf(tmp, 255, "CELL %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n",
+                         a*1e10, b*1e10, c*1e10,
+                         rad2deg(alpha), rad2deg(beta), rad2deg(gamma));
+		mosflm_sendline(tmp, mosflm);
+		break;
+
+	case 3 :
+		sg = cell_get_spacegroup(mosflm->target_cell);
+		/* Remove white space from space group */
+		j = 0;
+		for ( i=0; i<strlen(sg); i++ ) {
+			if (sg[i] != ' ') {
+				symm[j++] = sg[i];
+			}
 		}
+		symm[j] = '\0';
+		snprintf(tmp, 255, "SYMM %s\n", symm);
+		mosflm_sendline(tmp, mosflm);
+		break;
+
+	case 4 :
+		mosflm_sendline("DISTANCE 67.8\n", mosflm);
+		break;
+
+	case 5 :
+		mosflm_sendline("BEAM 0.0 0.0\n", mosflm);
+		break;
+
+	case 6 :
+		wavelength = image->lambda*1e10;
+		snprintf(tmp, 255, "WAVELENGTH %10.5f\n", wavelength);
+		mosflm_sendline(tmp, mosflm);
+		break;
+
+	case 7 :
+		snprintf(tmp, 255, "NEWMAT %s\n", mosflm->newmatfile);
+		mosflm_sendline(tmp, mosflm);
+		break;
+
+	case 8 :
+		snprintf(tmp, 255, "IMAGE %s phi 0 0\n", mosflm->imagefile);
+		mosflm_sendline(tmp, mosflm);
+		break;
+
+	case 9 :
+		snprintf(tmp, 255, "AUTOINDEX DPS FILE %s IMAGE 1\n",
+		         mosflm->sptfile);
+		mosflm_sendline(tmp, mosflm);
+		break;
+
+	case 10 :
+		mosflm_sendline("GO\n", mosflm);
+		break;
+
+	default:
+		mosflm_sendline("exit\n", mosflm);
+		return;
+
 	}
-	symm[j] = '\0';
-
-
-	/* build a script to run mosflm */
-	sprintf(mos_cmd,"%s","ipmosflm << eof-mosflm > /dev/null\n");
-	sprintf(mos_cmd,"%s%s",mos_cmd,
-	                     "DETECTOR ROTATION HORIZONTAL ANTICLOCKWISE"
-	                     " ORIGIN LL FAST HORIZONTAL RECTANGULAR\n");
-	sprintf(mos_cmd,"%sCELL %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n",
-	                   mos_cmd,
-	                   a*1e10,b*1e10,c*1e10,
-	                   rad2deg(alpha),rad2deg(beta),rad2deg(gamma));
-	sprintf(mos_cmd,"%sSYMM %s\n",mos_cmd,symm);
-	sprintf(mos_cmd,"%sDISTANCE %8.4f\n",mos_cmd,67.8);
-	sprintf(mos_cmd,"%sBEAM %8.4f %8.4f\n",mos_cmd,0.0,0.0);
-	sprintf(mos_cmd,"%sWAVELENGTH %10.5f\n",mos_cmd,wavelength);
-	sprintf(mos_cmd,"%sNEWMAT %s\n",mos_cmd,newmatfile);
-	sprintf(mos_cmd,"%sIMAGE xfel-%i_001.img phi 0 0\n",mos_cmd,image->id);
-	sprintf(mos_cmd,"%sAUTOINDEX DPS FILE xfel-%i.spt IMAGE 1\n",
-	                      mos_cmd,image->id);
-	sprintf(mos_cmd,"%sGO\n",mos_cmd);
-	sprintf(mos_cmd,"%s%s",mos_cmd,"eof-mosflm\n");
-
-	/* remove the previous NEWMAT file prior to running mosflm */
-	remove(newmatfile);
-
-	/* Run the mosflm script */
-	pid = fork();
-	if ( !( (pid != 0) && (pid != -1) ) ) {
-		if ( pid == -1 ) {
-			ERROR("fork() failed.\n");
-		} else {
 
-			/* Forked successfully, child process */
-			if ( system(mos_cmd) ) {
-				ERROR("MOSFLM execution failed.\n");
-				exit(0);
-				return;
+	mosflm->step++;
+}
+
+
+static int mosflm_readable(struct image *image, struct mosflm_data *mosflm)
+{
+	int rval;
+	int no_string = 0;
+
+	rval = read(mosflm->pty, mosflm->rbuffer+mosflm->rbufpos,
+	            mosflm->rbuflen-mosflm->rbufpos);
+
+	if ( (rval == -1) || (rval == 0) ) return 1;
+
+	mosflm->rbufpos += rval;
+	assert(mosflm->rbufpos <= mosflm->rbuflen);
+
+	while ( (!no_string) && (mosflm->rbufpos > 0) ) {
+
+		int i;
+		int block_ready = 0;
+		MOSFLMInputType type = MOSFLM_INPUT_NONE;
+
+		/* See if there's a full line in the buffer yet */
+		for ( i=0; i<mosflm->rbufpos-1; i++ ) {
+			/* Means the last value looked at is rbufpos-2 */
+
+			/* Is there a prompt in the buffer? */
+			if ( (i+7 <= mosflm->rbufpos)
+			  && (!strncmp(mosflm->rbuffer+i, "MOSFLM => ", 10)) ) {
+				block_ready = 1;
+				type = MOSFLM_INPUT_PROMPT;
+				break;
+			}
+
+			if ( (mosflm->rbuffer[i] == '\r')
+			  && (mosflm->rbuffer[i+1] == '\n') ) {
+				block_ready = 1;
+				type = MOSFLM_INPUT_LINE;
+				break;
+			}
+
+		}
+
+		if ( block_ready ) {
+
+			unsigned int new_rbuflen;
+			unsigned int endbit_length;
+			char *block_buffer = NULL;
+
+			switch ( type ) {
+
+			case MOSFLM_INPUT_LINE :
+
+				block_buffer = malloc(i+1);
+				memcpy(block_buffer, mosflm->rbuffer, i);
+				block_buffer[i] = '\0';
+
+				if ( block_buffer[0] == '\r' ) {
+					memmove(block_buffer, block_buffer+1, i);
+				}
+
+				mosflm_parseline(block_buffer, image, mosflm);
+				free(block_buffer);
+				endbit_length = i+2;
+				break;
+
+			case MOSFLM_INPUT_PROMPT :
+
+				mosflm_send_next(image, mosflm);
+				endbit_length = i+7;
+				break;
+
+			default :
+
+				/* Obviously, this never happens :) */
+				ERROR("Unrecognised MOSFLM input mode! "
+				      "I don't know how to understand MOSFLM\n");
+				return 1;
+
 			}
 
-			exit(0);
+			/* Now the block's been parsed, it should be
+			 * forgotten about */
+			memmove(mosflm->rbuffer,
+			        mosflm->rbuffer + endbit_length,
+			        mosflm->rbuflen - endbit_length);
+
+			/* Subtract the number of bytes removed */
+			mosflm->rbufpos = mosflm->rbufpos
+			                       - endbit_length;
+			new_rbuflen = mosflm->rbuflen - endbit_length;
+			if ( new_rbuflen == 0 ) new_rbuflen = 256;
+			mosflm->rbuffer = realloc(mosflm->rbuffer,
+			                               new_rbuflen);
+			mosflm->rbuflen = new_rbuflen;
+
+		} else {
+
+			if ( mosflm->rbufpos==mosflm->rbuflen ) {
+
+				/* More buffer space is needed */
+				mosflm->rbuffer = realloc(
+				                    mosflm->rbuffer,
+				                    mosflm->rbuflen + 256);
+				mosflm->rbuflen = mosflm->rbuflen
+				                                          + 256;
+				/* The new space gets used at the next
+				 * read, shortly... */
+
+			}
+			no_string = 1;
 
 		}
+
+	}
+
+	return 0;
+}
+
+
+void run_mosflm(struct image *image, UnitCell *cell)
+{
+	int fail;
+	struct mosflm_data *mosflm;
+	unsigned int opts;
+	int status;
+	int rval;
+
+	mosflm = malloc(sizeof(struct mosflm_data));
+	if ( mosflm == NULL ) {
+		ERROR("Couldn't allocate memory for MOSFLM data.\n");
+		return;
 	}
-	waitpid(pid, &r, 0);
 
-	/* Read the mosflm NEWMAT file and set cell candidate */
-	/* Existence of this file means possible success. Pretty shady. */
-	fail = read_newmat(newmatfile,image);
-	if (fail) {
-		printf("Failed to read mosflm NEWMAT file.\n");
+	mosflm->target_cell = cell;
+
+	snprintf(mosflm->imagefile, 127, "xfel-%i_001.img", image->id);
+	write_img(image, mosflm->imagefile); /* Dummy image */
+
+	snprintf(mosflm->sptfile, 127, "xfel-%i_001.spt", image->id);
+	write_spt(image, mosflm->sptfile);
+
+	snprintf(mosflm->newmatfile, 127, "xfel-%i.newmat", image->id);
+
+	mosflm->pid = forkpty(&mosflm->pty, NULL, NULL, NULL);
+	if ( mosflm->pid == -1 ) {
+		ERROR("Failed to fork for MOSFLM\n");
 		return;
 	}
+	if ( mosflm->pid == 0 ) {
+
+		/* Child process: invoke MOSFLM */
+		struct termios t;
+
+		/* Turn echo off */
+		tcgetattr(STDIN_FILENO, &t);
+		t.c_lflag &= ~(ECHO | ECHOE | ECHOK | ECHONL);
+		tcsetattr(STDIN_FILENO, TCSANOW, &t);
+
+		execlp("mosflm", "", (char *)NULL);
+		ERROR("Failed to invoke MOSFLM.\n");
+		_exit(0);
+
+	}
+
+	mosflm->rbuffer = malloc(256);
+	mosflm->rbuflen = 256;
+	mosflm->rbufpos = 0;
+
+	/* Set non-blocking */
+	opts = fcntl(mosflm->pty, F_GETFL);
+	fcntl(mosflm->pty, F_SETFL, opts | O_NONBLOCK);
+
+	mosflm->step = 1;	/* This starts the "initialisation" procedure */
+
+	do {
+
+		fd_set fds;
+		struct timeval tv;
+		int sval;
+
+		FD_ZERO(&fds);
+		FD_SET(mosflm->pty, &fds);
+
+		tv.tv_sec = 10;
+		tv.tv_usec = 0;
+
+		sval = select(mosflm->pty+1, &fds, NULL, NULL, &tv);
+
+		if ( sval == -1 ) {
+			int err = errno;
+			ERROR("select() failed: %s\n", strerror(err));
+			rval = 1;
+		} else if ( sval != 0 ) {
+			rval = mosflm_readable(image, mosflm);
+		} else {
+			ERROR("No response from MOSFLM..\n");
+			rval = 1;
+		}
+
+	} while ( !rval );
+
+	close(mosflm->pty);
+	free(mosflm->rbuffer);
+	wait(&status);
+
+	/* Read the mosflm NEWMAT file and set cell candidate *
+	 * Existence of this file means possible success. Pretty shady. */
+	fail = read_newmat(mosflm->newmatfile, image);
+	if ( fail ) {
+		ERROR("Failed to read MOSFLM's NEWMAT file.\n");
+		return;
+	}
+
+	free(mosflm);
 }