interpolate and various sizes build

diff --git a/.bzrignore b/.bzrignore
index ac0ffb3eec608d0d84143ae42005b5643dd9deed..2f356c250c36c1344c08d41dd7c870ffac038710 100644 (file)
--- 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 9c82be7b10a0e3dc4eff139e1a01b0ade6d3fe8f..42ebbd59caca045940666ede5f5396a780292cf5 100644 (file)
--- 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 34381e890494a489c95e0161856dccdc598a6ce6..62b981efa930fd19bbb3eab3bb94dd0ba0d14d69 100644 (file)
--- 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 7e0e86d7622da0698d992bdf94af0dfef3e4676e..99793574edf25789f86f9b092ac1d1f4682e4fb7 100644 (file)
--- 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 (file)
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 cc86803594841154686a48adb5c77aff1bb5b316..0121cad460d6ba29972c59471e44452b1527a756 100644 (file)
--- 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"