From: Ian Jackson <ian@davenant.relativity.greenend.org.uk>
Date: Fri, 18 Jan 2008 20:33:56 +0000 (+0000)
Subject: interpolate and various sizes build
X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?p=moebius2.git;a=commitdiff_plain;h=97ef96dd318b684d829ac950403582788df427d8

interpolate and various sizes build
---

diff --git a/.bzrignore b/.bzrignore
index ac0ffb3..2f356c2 100644
--- a/.bzrignore
+++ b/.bzrignore
@@ -1,6 +1,7 @@
 minimise
 primer
-view+[1-9][1-9]
+view-[1-9][1-9]
+interpolate-[1-9][1-9]
 *.d
 *.tmp
 *.new
diff --git a/Makefile b/Makefile
index 9c82be7..42ebbd5 100644
--- a/Makefile
+++ b/Makefile
@@ -2,7 +2,8 @@
 
 VIEWDIMS=33 44 55
 TARGETS= minimise primer lumpy.cfm sgtatham.cfm ring.cfm \
-	$(addprefix view+, $(VIEWDIMS))
+	interpolate-44 \
+	$(addprefix view-, $(VIEWDIMS))
 SGTATHAM=sgtatham
 
 CWARNS=	-Wall -Wwrite-strings -Wpointer-arith -Werror -Wshadow
@@ -37,20 +38,24 @@ lumpy.cfm: oldmoebius-converter prime.data ../moebius/ins-new ../moebius/a.out
 ring.cfm: oldmoebius-converter prime.data /dev/null ../moebius/a.out
 		./$^ -o$@
 
+view-%:		view+%.o mgraph+%.o common.o
+		$(CC) $(CFLAGS) -o $@ $^ $(LIBGSL) -L/usr/X11R6/lib -lX11
 
-# this ridiculous repetition is due to make being too lame
+interpolate-%:	interpolate+%.o mgraph+%.o common.o
+		$(CC) $(CFLAGS) -o $@ $^ $(LIBGSL)
 
-view+%:		view+%.o mgraph+%.o common.o
-		$(CC) $(CFLAGS) -o $@ $^ $(LIBGSL) -L/usr/X11R6/lib -lX11
+# this ridiculous repetition is due to make being too lame
 
-view+%.o:	view.c
+view+%.o: view.c
 		$(CC) -c $(CPPFLAGS) $(CFLAGS) -DDEFBITS=$* $< -o $@
 
-mgraph+%.o:	mgraph.c
+mgraph+%.o: mgraph.c
 		$(CC) -c $(CPPFLAGS) $(CFLAGS) -DDEFBITS=$* $< -o $@
 
-.PRECIOUS: view+%.o mgraph+%.o
+interpolate+%.o: interpolate.c
+		$(CC) -c $(CPPFLAGS) $(CFLAGS) -DDEFBITS=$* $< -o $@
 
+.PRECIOUS: view+%.o mgraph+%.o interpolate+%.o
 
 clean:
 		rm -f prime.data $(TARGETS)
@@ -60,4 +65,6 @@ clean:
 realclean:	clean
 		rm -f best
 
+%.d:
+
 -include *.d
diff --git a/common.c b/common.c
index 34381e8..62b981e 100644
--- a/common.c
+++ b/common.c
@@ -61,4 +61,5 @@ void gsldie(int l, const char *what, int status) {
 }
 
 void diee(const char *what) { perror(what); exit(16); }
+void fail(const char *emsg) { fputs(emsg,stderr); exit(12); }
 void flushoutput(void) { if (fflush(stdout)||ferror(stdout)) diee("stdout"); }
diff --git a/common.h b/common.h
index 7e0e86d..9979357 100644
--- a/common.h
+++ b/common.h
@@ -9,6 +9,9 @@
 #define _GNU_SOURCE
 #endif
 
+#include <sys/types.h>
+#include <sys/stat.h>
+
 #include <math.h>
 #include <float.h>
 #include <limits.h>
@@ -23,6 +26,7 @@
 #include <gsl/gsl_matrix.h>
 #include <gsl/gsl_blas.h>
 #include <gsl/gsl_linalg.h>
+#include <gsl/gsl_interp.h>
 
 #define D3 3
 
@@ -36,6 +40,7 @@ double dotprod(const double a[D3], const double b[D3]);
 
 void flushoutput(void);
 void diee(const char *what);
+void fail(const char *emsg);
 
 void libdie(const char *lib, int l, const char *str);
 #define XA(w) ((w) ? (void)0 : libdie("X", __LINE__, #w))
diff --git a/interpolate.c b/interpolate.c
new file mode 100644
index 0000000..6f519e4
--- /dev/null
+++ b/interpolate.c
@@ -0,0 +1,136 @@
+/*
+ * increases the resolution of a conformation by interpolation
+ */
+
+#include "mgraph.h"
+
+/* filled in by characterise_input: */
+static int shift, oldxbits, oldybits, oldx, oldy, oldsz, inc;
+
+/* filled in by read_input: */
+static struct Vertices all;
+
+static void characterise_input(void) {
+  struct stat stab;
+  int r;
+  
+  r= fstat(0,&stab);  if (r) diee("fstat input to find length");
+
+  if (!stab.st_size || stab.st_size % sizeof(double) ||
+      stab.st_size > INT_MAX)
+    fail("input file is not reasonable whole number of doubles\n");
+
+  oldsz= stab.st_size / sizeof(double);
+  for (shift=1;
+       shift > XBITS+1 && shift > YBITS+1;
+       shift++) {
+    oldxbits= XBITS-1;
+    oldybits= YBITS-1;
+    oldx=  1<<oldxbits;
+    oldy= (1<<oldybits)-1;
+    if (oldx*oldy == oldsz) goto found;
+  }
+  fail("input file size cannot be interpolated to target file size\n");
+
+ found:
+  inc= 1<<shift;
+}
+
+static void read_input(void) {
+  int x,y;
+  
+  for (y=0; y<Y; y+=inc)
+    for (x=0; x<X; x+=inc) {
+      errno= 0;
+      if (fread(all.a[(y << YSHIFT) | x], sizeof(double), D3, stdin) != D3)
+	diee("read input");
+    }
+}
+
+  /* We use GSL's interpolation functions.  Each xa array is simple
+   * integers, and we do the interpolation three times for each line
+   * of topologically colinear vertices (once, separately, for each of
+   * the three spatial coordinates in the ya array).  The `horizontal'
+   * lines are done with periodic boundary conditions, obviously.
+   *
+   * This directly gives the required additional vertices:
+   *
+   *         ,O.                   O  original vertex
+   *        /   \                  *  new vertex
+   *      ,* - - *.                dashed lines are new edges,
+   *     /  `   '  \                not directly represented
+   *    O----`*'----O               or manipulated by this program
+   *
+   */
+
+#define NEXTV (v= EDGE_END2(v,direction))
+
+static void interpolate_line(int startvertex,
+			     int direction /* edge number */,
+			     int norig) {
+  double xa[norig+1], ya[D3][norig+1], n;
+  const gsl_interp_type *it;
+  gsl_interp *interp;
+  gsl_interp_accel *accel;
+  int i, v, k;
+
+  for (i=0, v=startvertex;
+       ;
+       i++, NEXTV, assert(v>=0), NEXTV) {
+    if (v<0)
+      break;
+    assert(i <= norig);
+    xa[i]= i;
+    K ya[k][i]= all.a[v][k];
+    if (v==startvertex)
+      break;
+  }
+  n= i;
+  it= v>=0 ? gsl_interp_cspline_periodic : gsl_interp_cspline;
+
+  K {
+    interp= gsl_interp_alloc(it,n); if (!interp) diee("gsl_interp_alloc");
+    accel= gsl_interp_accel_alloc(); if(!accel) diee("gsl_interp_accel_alloc");
+    GA( gsl_interp_init(interp,xa,ya[k],n) );
+
+    for (i=0, v=startvertex;
+	 i<n-1;
+	 i++, NEXTV) {
+      assert(v>=0); NEXTV; assert(v>=0);
+      GA( gsl_interp_eval_e(interp,xa,ya[k], i+0.5, accel, &all.a[v][k]) );
+    }
+    
+    gsl_interp_accel_free(accel);
+    gsl_interp_free(interp);
+  }
+}
+	
+static void interpolate(void) {
+  int x,y;
+  
+  if (shift!=1) fail("only interpolation by factor of 2 supported\n");
+
+  for (y=0; y<Y; y+=inc) {
+    interpolate_line(y<<YSHIFT, 0, oldx);
+  }
+  for (x=0; x<X; x+=inc) {
+    interpolate_line(x, 5, oldy);
+    interpolate_line(x, 4, oldy);
+  }
+}
+
+static void write_output(void) {
+  if (fwrite(all.a,sizeof(all.a),1,stdout) != 1) diee("fwrite");
+}
+
+int main(int argc, const char **argv) {
+  if (argc!=1 || isatty(0) || isatty(1))
+    { fputs("usage: interpolate <INPUT.cfm >OUTPUT.cfm\n",stderr); exit(8); }
+
+  characterise_input();
+  read_input();
+  interpolate();
+  write_output();
+  if (fclose(stdout)) diee("fclose stdout");
+  return 0;
+}
diff --git a/view.c b/view.c
index cc86803..0121cad 100644
--- a/view.c
+++ b/view.c
@@ -5,8 +5,6 @@
 #include <X11/Xlib.h>
 #include <X11/Xutil.h>
 
-#include <sys/types.h>
-#include <sys/stat.h>
 #include <sys/poll.h>
 
 #include "mgraph.h"