UnitCell: Store all representations once they're calculated

Previously, the "getter" functions would re-calculate the requested
representation every time they were called.  This could mean doing a
matrix inversion in the middle of a tight loop, wasting loads of time.
Now, it stores the calculated values and returns them directly next
time.  Setting the parameters invalidates the values for all
representations other than the one given.

The cost of doing this is that the cell can no longer be "const" in the
getter functions.  This tracked through some other code, but nothing too
severe.

6 files changed, 174 insertions, 195 deletions

diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index 93dea212..6b5d7717 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -333,8 +333,6 @@ void cell_print_full(UnitCell *cell)
 		       rad2deg(angle_between(asx, asy, asz, csx, csy, csz)),
 		       rad2deg(angle_between(asx, asy, asz, bsx, bsy, bsz)));
 
-		STATUS("Cell representation is %s.\n", cell_rep(cell));
-
 	}
 }
 
@@ -910,7 +908,7 @@ static int get_angle_rad(char **bits, int nbits, double *pl)
  * Writes \p cell to \p fh, in CrystFEL unit cell file format
  *
  */
-void write_cell(const UnitCell *cell, FILE *fh)
+void write_cell(UnitCell *cell, FILE *fh)
 {
 	double a, b, c, al, be, ga;
 	LatticeType lt;
diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h
index 0a110bf2..5789d606 100644
--- a/libcrystfel/src/cell-utils.h
+++ b/libcrystfel/src/cell-utils.h
@@ -58,7 +58,7 @@ extern void cell_print_full(UnitCell *cell);
 
 extern UnitCell *load_cell_from_pdb(const char *filename);
 extern UnitCell *load_cell_from_file(const char *filename);
-extern void write_cell(const UnitCell *cell, FILE *fh);
+extern void write_cell(UnitCell *cell, FILE *fh);
 
 extern int cell_is_sensible(UnitCell *cell);
 
diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c
index 23a5e46b..50a28470 100644
--- a/libcrystfel/src/cell.c
+++ b/libcrystfel/src/cell.c
@@ -35,6 +35,7 @@
 #endif
 
 #include <math.h>
+#include <assert.h>
 #include <stdlib.h>
 #include <stdio.h>
 #include <string.h>
@@ -54,19 +55,14 @@
  */
 
 
-typedef enum {
-	CELL_REP_CRYST,
-	CELL_REP_CART,
-	CELL_REP_RECIP
-} CellRepresentation;
-
 struct _unitcell {
 
-	CellRepresentation rep;
-
-	int have_parameters;
+	LatticeType  lattice_type;
+	char         centering;
+	char         unique_axis;
 
 	/* Crystallographic representation */
+	int have_cryst;
 	double a;	/* m */
 	double b;	/* m */
 	double c;	/* m */
@@ -75,18 +71,16 @@ struct _unitcell {
 	double gamma;	/* Radians */
 
 	/* Cartesian representation */
+	int have_cart;
 	double ax;	double bx;	double cx;
 	double ay;	double by;	double cy;
 	double az;	double bz;	double cz;
 
 	/* Cartesian representation of reciprocal axes */
+	int have_recip;
 	double axs;	double bxs;	double cxs;
 	double ays;	double bys;	double cys;
 	double azs;	double bzs;	double czs;
-
-	LatticeType  lattice_type;
-	char         centering;
-	char         unique_axis;
 };
 
 typedef enum {
@@ -127,12 +121,13 @@ UnitCell *cell_new()
 	cell->beta = 0.0;
 	cell->gamma = 0.0;
 
-	cell->rep = CELL_REP_CRYST;
+	cell->have_cryst = 0;
+	cell->have_cart = 0;
+	cell->have_recip = 0;
 
 	cell->lattice_type = L_TRICLINIC;
 	cell->centering = 'P';
 	cell->unique_axis = '?';
-	cell->have_parameters = 0;
 
 	return cell;
 }
@@ -160,7 +155,9 @@ void cell_free(UnitCell *cell)
 int cell_has_parameters(const UnitCell *cell)
 {
 	if ( cell == NULL ) return 0;
-	return cell->have_parameters;
+	return (cell->have_cryst > 0)
+	    || (cell->have_cart > 0 )
+	    || (cell->have_recip > 0);
 }
 
 
@@ -176,8 +173,9 @@ void cell_set_parameters(UnitCell *cell, double a, double b, double c,
 	cell->beta = beta;
 	cell->gamma = gamma;
 
-	cell->rep = CELL_REP_CRYST;
-	cell->have_parameters = 1;
+	cell->have_cryst = 1;
+	cell->have_cart = 0;
+	cell->have_recip = 0;
 }
 
 
@@ -192,8 +190,9 @@ void cell_set_cartesian(UnitCell *cell,
 	cell->bx = bx;  cell->by = by;  cell->bz = bz;
 	cell->cx = cx;  cell->cy = cy;  cell->cz = cz;
 
-	cell->rep = CELL_REP_CART;
-	cell->have_parameters = 1;
+	cell->have_cryst = 0;
+	cell->have_cart = 1;
+	cell->have_recip = 0;
 }
 
 
@@ -223,8 +222,9 @@ UnitCell *cell_new_from_reciprocal_axes(struct rvec as, struct rvec bs,
 	cell->bxs = bs.u;  cell->bys = bs.v;  cell->bzs = bs.w;
 	cell->cxs = cs.u;  cell->cys = cs.v;  cell->czs = cs.w;
 
-	cell->rep = CELL_REP_RECIP;
-	cell->have_parameters = 1;
+	cell->have_cryst = 0;
+	cell->have_cart = 0;
+	cell->have_recip = 1;
 
 	return cell;
 }
@@ -241,8 +241,9 @@ UnitCell *cell_new_from_direct_axes(struct rvec a, struct rvec b, struct rvec c)
 	cell->bx = b.u;  cell->by = b.v;  cell->bz = b.w;
 	cell->cx = c.u;  cell->cy = c.v;  cell->cz = c.w;
 
-	cell->rep = CELL_REP_CART;
-	cell->have_parameters = 1;
+	cell->have_cryst = 0;
+	cell->have_cart = 1;
+	cell->have_recip = 0;
 
 	return cell;
 }
@@ -268,8 +269,9 @@ void cell_set_reciprocal(UnitCell *cell,
 	cell->bxs = bsx;  cell->bys = bsy;  cell->bzs = bsz;
 	cell->cxs = csx;  cell->cys = csy;  cell->czs = csz;
 
-	cell->rep = CELL_REP_RECIP;
-	cell->have_parameters = 1;
+	cell->have_cryst = 0;
+	cell->have_cart = 0;
+	cell->have_recip = 1;
 }
 
 
@@ -291,31 +293,25 @@ void cell_set_unique_axis(UnitCell *cell, char unique_axis)
 }
 
 
-/************************* Getter helper functions ****************************/
+/************************* Conversion functions ****************************/
 
-static int cell_crystallographic_to_cartesian(const UnitCell *cell,
-                                              double *ax, double *ay, double *az,
-                                              double *bx, double *by, double *bz,
-                                              double *cx, double *cy, double *cz)
+static void crystallographic_to_cartesian(UnitCell *cell)
 {
 	double tmp, V, cosalphastar, cstar;
 
-	if ( !cell->have_parameters ) {
-		ERROR("Unit cell has unspecified parameters.\n");
-		return 1;
-	}
+	assert(cell->have_cryst == 1);
 
 	/* Firstly: Get a in terms of x, y and z
 	 * +a (cryst) is defined to lie along +x (cart) */
-	*ax = cell->a;
-	*ay = 0.0;
-	*az = 0.0;
+	cell->ax = cell->a;
+	cell->ay = 0.0;
+	cell->az = 0.0;
 
 	/* b in terms of x, y and z
 	 * b (cryst) is defined to lie in the xy (cart) plane */
-	*bx = cell->b*cos(cell->gamma);
-	*by = cell->b*sin(cell->gamma);
-	*bz = 0.0;
+	cell->bx = cell->b*cos(cell->gamma);
+	cell->by = cell->b*sin(cell->gamma);
+	cell->bz = 0.0;
 
 	tmp = cos(cell->alpha)*cos(cell->alpha)
 		+ cos(cell->beta)*cos(cell->beta)
@@ -329,21 +325,41 @@ static int cell_crystallographic_to_cartesian(const UnitCell *cell,
 	cstar = (cell->a * cell->b * sin(cell->gamma))/V;
 
 	/* c in terms of x, y and z */
-	*cx = cell->c*cos(cell->beta);
-	*cy = -cell->c*sin(cell->beta)*cosalphastar;
-	*cz = 1.0/cstar;
+	cell->cx = cell->c*cos(cell->beta);
+	cell->cy = -cell->c*sin(cell->beta)*cosalphastar;
+	cell->cz = 1.0/cstar;
 
-	return 0;
+	cell->have_cart = 1;
+}
+
+
+static void cartesian_to_crystallographic(UnitCell *cell)
+{
+	assert(cell->have_cart);
+
+	/* Convert cartesian -> crystallographic */
+	cell->a = modulus(cell->ax, cell->ay, cell->az);
+	cell->b = modulus(cell->bx, cell->by, cell->bz);
+	cell->c = modulus(cell->cx, cell->cy, cell->cz);
+
+	cell->alpha = angle_between(cell->bx, cell->by, cell->bz,
+	                            cell->cx, cell->cy, cell->cz);
+	cell->beta = angle_between(cell->ax, cell->ay, cell->az,
+	                           cell->cx, cell->cy, cell->cz);
+	cell->gamma = angle_between(cell->ax, cell->ay, cell->az,
+	                            cell->bx, cell->by, cell->bz);
+
+	cell->have_cryst = 1;
 }
 
 
 /* Why yes, I do enjoy long argument lists...! */
-static int cell_invert(double ax, double ay, double az,
-                       double bx, double by, double bz,
-                       double cx, double cy, double cz,
-                       double *asx, double *asy, double *asz,
-                       double *bsx, double *bsy, double *bsz,
-                       double *csx, double *csy, double *csz)
+static int invert(double ax, double ay, double az,
+                  double bx, double by, double bz,
+                  double cx, double cy, double cz,
+                  double *asx, double *asy, double *asz,
+                  double *bsx, double *bsy, double *bsz,
+                  double *csx, double *csy, double *csz)
 {
 	int s;
 	gsl_matrix *m;
@@ -413,157 +429,143 @@ static int cell_invert(double ax, double ay, double az,
 }
 
 
+static int reciprocal_to_cartesian(UnitCell *cell)
+{
+	assert(cell->have_recip);
+
+	if ( invert(cell->axs, cell->ays, cell->azs,
+	            cell->bxs, cell->bys, cell->bzs,
+	            cell->cxs, cell->cys, cell->czs,
+	            &cell->ax, &cell->ay, &cell->az,
+	            &cell->bx, &cell->by, &cell->bz,
+	            &cell->cx, &cell->cy, &cell->cz) ) return 1;
+
+	cell->have_cart = 1;
+	return 0;
+}
+
+
+static int cartesian_to_reciprocal(UnitCell *cell)
+{
+	assert(cell->have_cart);
+
+	if ( invert(cell->ax, cell->ay, cell->az,
+	            cell->bx, cell->by, cell->bz,
+	            cell->cx, cell->cy, cell->cz,
+	            &cell->axs, &cell->ays, &cell->azs,
+	            &cell->bxs, &cell->bys, &cell->bzs,
+	            &cell->cxs, &cell->cys, &cell->czs) ) return 1;
+
+	cell->have_recip = 1;
+	return 0;
+}
+
+
 /********************************** Getters ***********************************/
 
-int cell_get_parameters(const UnitCell *cell, double *a, double *b, double *c,
+int cell_get_parameters(UnitCell *cell,
+                        double *a, double *b, double *c,
                         double *alpha, double *beta, double *gamma)
 {
-	double ax, ay, az, bx, by, bz, cx, cy, cz;
-
 	if ( cell == NULL ) return 1;
 
-	if ( !cell->have_parameters ) {
+	if ( cell->have_cryst ) {
+
+		/* Nothing to do */
+
+	} else if ( cell->have_cart ) {
+
+		cartesian_to_crystallographic(cell);
+
+	} else if ( cell->have_recip ) {
+
+		if ( reciprocal_to_cartesian(cell) ) return 1;
+		cartesian_to_crystallographic(cell);
+
+	} else {
+
 		ERROR("Unit cell has unspecified parameters.\n");
 		return 1;
-	}
 
-	switch ( cell->rep ) {
-
-		case CELL_REP_CRYST:
-		/* Direct response */
-		*a = cell->a;
-		*b = cell->b;
-		*c = cell->c;
-		*alpha = cell->alpha;
-		*beta = cell->beta;
-		*gamma = cell->gamma;
-		return 0;
-
-		case CELL_REP_CART:
-		/* Convert cartesian -> crystallographic */
-		*a = modulus(cell->ax, cell->ay, cell->az);
-		*b = modulus(cell->bx, cell->by, cell->bz);
-		*c = modulus(cell->cx, cell->cy, cell->cz);
-
-		*alpha = angle_between(cell->bx, cell->by, cell->bz,
-		                       cell->cx, cell->cy, cell->cz);
-		*beta = angle_between(cell->ax, cell->ay, cell->az,
-		                      cell->cx, cell->cy, cell->cz);
-		*gamma = angle_between(cell->ax, cell->ay, cell->az,
-		                       cell->bx, cell->by, cell->bz);
-		return 0;
-
-		case CELL_REP_RECIP:
-		/* Convert reciprocal -> crystallographic.
-                 * Start by converting reciprocal -> cartesian */
-		if ( cell_invert(cell->axs, cell->ays, cell->azs,
-		                 cell->bxs, cell->bys, cell->bzs,
-		                         cell->cxs, cell->cys, cell->czs,
-		                         &ax, &ay, &az,
-		                         &bx, &by, &bz,
-		                         &cx, &cy, &cz) ) return 1;
-
-		/* Now convert cartesian -> crystallographic */
-		*a = modulus(ax, ay, az);
-		*b = modulus(bx, by, bz);
-		*c = modulus(cx, cy, cz);
-
-		*alpha = angle_between(bx, by, bz, cx, cy, cz);
-		*beta = angle_between(ax, ay, az, cx, cy, cz);
-		*gamma = angle_between(ax, ay, az, bx, by, bz);
-		return 0;
 	}
 
-	return 1;
+	*a = cell->a;
+	*b = cell->b;
+	*c = cell->c;
+	*alpha = cell->alpha;
+	*beta = cell->beta;
+	*gamma = cell->gamma;
+	return 0;
 }
 
 
-int cell_get_cartesian(const UnitCell *cell,
+int cell_get_cartesian(UnitCell *cell,
                        double *ax, double *ay, double *az,
                        double *bx, double *by, double *bz,
                        double *cx, double *cy, double *cz)
 {
 	if ( cell == NULL ) return 1;
 
-	if ( !cell->have_parameters ) {
+	if ( cell->have_cart ) {
+
+		/* Nothing to do */
+
+	} else if ( cell->have_recip ) {
+
+		/* NB recip->cart has priority over
+		 * cryst->cart, to preserve orientation */
+		reciprocal_to_cartesian(cell);
+
+	} else if ( cell->have_cryst ) {
+
+		crystallographic_to_cartesian(cell);
+
+	} else {
+
 		ERROR("Unit cell has unspecified parameters.\n");
 		return 1;
-	}
-
-	switch ( cell->rep ) {
-
-		case CELL_REP_CRYST:
-		/* Convert crystallographic -> cartesian. */
-		return cell_crystallographic_to_cartesian(cell,
-		                                          ax, ay, az,
-		                                          bx, by, bz,
-		                                          cx, cy, cz);
-
-		case CELL_REP_CART:
-		/* Direct response */
-		*ax = cell->ax;  *ay = cell->ay;  *az = cell->az;
-		*bx = cell->bx;  *by = cell->by;  *bz = cell->bz;
-		*cx = cell->cx;  *cy = cell->cy;  *cz = cell->cz;
-		return 0;
-
-		case CELL_REP_RECIP:
-		/* Convert reciprocal -> cartesian */
-		return cell_invert(cell->axs, cell->ays, cell->azs,
-		                   cell->bxs, cell->bys, cell->bzs,
-		                   cell->cxs, cell->cys, cell->czs,
-		                   ax, ay, az, bx, by, bz, cx, cy, cz);
 
 	}
 
-	return 1;
+	*ax = cell->ax;  *ay = cell->ay;  *az = cell->az;
+	*bx = cell->bx;  *by = cell->by;  *bz = cell->bz;
+	*cx = cell->cx;  *cy = cell->cy;  *cz = cell->cz;
+	return 0;
 }
 
 
-int cell_get_reciprocal(const UnitCell *cell,
+int cell_get_reciprocal(UnitCell *cell,
                         double *asx, double *asy, double *asz,
                         double *bsx, double *bsy, double *bsz,
                         double *csx, double *csy, double *csz)
 {
-	int r;
-	double ax, ay, az, bx, by, bz, cx, cy, cz;
-
 	if ( cell == NULL ) return 1;
 
-	if ( !cell->have_parameters ) {
+	if ( cell->have_recip ) {
+
+		/* Nothing to do */
+
+	} else if ( cell->have_cart ) {
+
+		/* NB cart->recip has priority over cryst->recip */
+		cartesian_to_reciprocal(cell);
+
+	} else if ( cell->have_cryst ) {
+
+		crystallographic_to_cartesian(cell);
+		cartesian_to_reciprocal(cell);
+
+	} else {
+
 		ERROR("Unit cell has unspecified parameters.\n");
 		return 1;
-	}
-
-	switch ( cell->rep ) {
-
-		case CELL_REP_CRYST:
-		/* Convert crystallographic -> reciprocal */
-		r = cell_crystallographic_to_cartesian(cell,
-		                                       &ax, &ay, &az,
-		                                       &bx, &by, &bz,
-		                                       &cx, &cy, &cz);
-		if ( r ) return r;
-		return cell_invert(ax, ay, az,bx, by, bz, cx, cy, cz,
-		                   asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
-
-		case CELL_REP_CART:
-		/* Convert cartesian -> reciprocal */
-		cell_invert(cell->ax, cell->ay, cell->az,
-		            cell->bx, cell->by, cell->bz,
-		            cell->cx, cell->cy, cell->cz,
-		            asx, asy, asz, bsx, bsy, bsz, csx, csy, csz);
-		return 0;
-
-		case CELL_REP_RECIP:
-		/* Direct response */
-		*asx = cell->axs;  *asy = cell->ays;  *asz = cell->azs;
-		*bsx = cell->bxs;  *bsy = cell->bys;  *bsz = cell->bzs;
-		*csx = cell->cxs;  *csy = cell->cys;  *csz = cell->czs;
-		return 0;
 
 	}
 
-	return 1;
+	*asx = cell->axs;  *asy = cell->ays;  *asz = cell->azs;
+	*bsx = cell->bxs;  *bsy = cell->bys;  *bsz = cell->bzs;
+	*csx = cell->cxs;  *csy = cell->cys;  *csz = cell->czs;
+	return 0;
 }
 
 
@@ -579,7 +581,7 @@ LatticeType cell_get_lattice_type(const UnitCell *cell)
 }
 
 
-struct g6 cell_get_G6(const UnitCell *cell)
+struct g6 cell_get_G6(UnitCell *cell)
 {
 	double a, b, c, al, be, ga;
 	struct g6 g;
@@ -600,25 +602,6 @@ char cell_get_unique_axis(const UnitCell *cell)
 }
 
 
-const char *cell_rep(UnitCell *cell)
-{
-	switch ( cell->rep ) {
-
-		case CELL_REP_CRYST:
-		return "crystallographic, direct space";
-
-		case CELL_REP_CART:
-		return "cartesian, direct space";
-
-		case CELL_REP_RECIP:
-		return "cartesian, reciprocal space";
-
-	}
-
-	return "unknown";
-}
-
-
 UnitCell *cell_transform_gsl_direct(UnitCell *in, gsl_matrix *m)
 {
 	gsl_matrix *c;
diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h
index bf5d689c..064dacd0 100644
--- a/libcrystfel/src/cell.h
+++ b/libcrystfel/src/cell.h
@@ -107,15 +107,15 @@ extern void cell_set_cartesian(UnitCell *cell,
 extern void cell_set_parameters(UnitCell *cell, double a, double b, double c,
 				double alpha, double beta, double gamma);
 
-extern int cell_get_parameters(const UnitCell *cell, double *a, double *b, double *c,
+extern int cell_get_parameters(UnitCell *cell, double *a, double *b, double *c,
                                double *alpha, double *beta, double *gamma);
 
-extern int cell_get_cartesian(const UnitCell *cell,
+extern int cell_get_cartesian(UnitCell *cell,
                               double *ax, double *ay, double *az,
                               double *bx, double *by, double *bz,
                               double *cx, double *cy, double *cz);
 
-extern int cell_get_reciprocal(const UnitCell *cell,
+extern int cell_get_reciprocal(UnitCell *cell,
                                double *asx, double *asy, double *asz,
                                double *bsx, double *bsy, double *bsz,
                                double *csx, double *csy, double *csz);
@@ -138,7 +138,7 @@ struct g6
 	double F;
 };
 
-extern struct g6 cell_get_G6(const UnitCell *cell);
+extern struct g6 cell_get_G6(UnitCell *cell);
 
 extern char cell_get_centering(const UnitCell *cell);
 extern void cell_set_centering(UnitCell *cell, char centering);
@@ -146,8 +146,6 @@ extern void cell_set_centering(UnitCell *cell, char centering);
 extern char cell_get_unique_axis(const UnitCell *cell);
 extern void cell_set_unique_axis(UnitCell *cell, char unique_axis);
 
-extern const char *cell_rep(UnitCell *cell);
-
 extern UnitCell *cell_transform_gsl_direct(UnitCell *in, gsl_matrix *m);
 
 extern UnitCell *cell_transform_rational(UnitCell *cell, RationalMatrix *m);
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
index 08f0e332..afaf40a5 100644
--- a/libcrystfel/src/stream.c
+++ b/libcrystfel/src/stream.c
@@ -1250,7 +1250,7 @@ Stream *stream_open_fd_for_write(int fd, const DataTemplate *dtempl)
 }
 
 
-void stream_write_target_cell(Stream *st, const UnitCell *cell)
+void stream_write_target_cell(Stream *st, UnitCell *cell)
 {
 	if ( cell == NULL ) return;
 	fprintf(st->fh, STREAM_CELL_START_MARKER"\n");
diff --git a/libcrystfel/src/stream.h b/libcrystfel/src/stream.h
index 3d57c7b6..f503f391 100644
--- a/libcrystfel/src/stream.h
+++ b/libcrystfel/src/stream.h
@@ -98,7 +98,7 @@ extern void stream_close(Stream *st);
 extern void stream_write_geometry_file(Stream *st,
                                        const char *geom_filename);
 extern void stream_write_target_cell(Stream *st,
-                                     const UnitCell *cell);
+                                     UnitCell *cell);
 extern void stream_write_commandline_args(Stream *st,
                                           int argc, char *argv[]);
 extern void stream_write_indexing_methods(Stream *st,