chiark / gitweb /
interpolates but not very well
authorIan Jackson <ian@davenant.greenend.org.uk>
Fri, 18 Jan 2008 23:22:44 +0000 (23:22 +0000)
committerIan Jackson <ian@davenant.greenend.org.uk>
Fri, 18 Jan 2008 23:22:44 +0000 (23:22 +0000)
.bzrignore
Makefile
interpolate.c
mgraph.c
mgraph.h

index 93c99e1..156d40a 100644 (file)
@@ -1,7 +1,7 @@
 minimise
 primer
-view-[1-9][1-9]
-interpolate-[1-9][1-9]
+view-[1-9]*
+interpolate-[1-9]*
 *.d
 *.tmp
 *.new
index a2409fc..dd4b218 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,8 +1,8 @@
 
 
-VIEWDIMS=33 44 55
+VIEWDIMS=33 64 125
 TARGETS= minimise primer lumpy.cfm sgtatham.cfm ring.cfm \
-       interpolate-44 \
+       interpolate-64 \
        $(addprefix view-, $(VIEWDIMS))
 SGTATHAM=sgtatham
 
@@ -52,13 +52,13 @@ interpolate-%:      interpolate+%.o mgraph+%.o common.o
 # this ridiculous repetition is due to make being too lame
 
 view+%.o: view.c
-               $(CC) -c $(CPPFLAGS) $(CFLAGS) -DDEFBITS=$* $< -o $@
+               $(CC) -c $(CPPFLAGS) $(CFLAGS) -DDEFSZ=$* $< -o $@
 
 mgraph+%.o: mgraph.c
-               $(CC) -c $(CPPFLAGS) $(CFLAGS) -DDEFBITS=$* $< -o $@
+               $(CC) -c $(CPPFLAGS) $(CFLAGS) -DDEFSZ=$* $< -o $@
 
 interpolate+%.o: interpolate.c
-               $(CC) -c $(CPPFLAGS) $(CFLAGS) -DDEFBITS=$* $< -o $@
+               $(CC) -c $(CPPFLAGS) $(CFLAGS) -DDEFSZ=$* $< -o $@
 
 .PRECIOUS: view+%.o mgraph+%.o interpolate+%.o
 
index d9d3910..1de2bea 100644 (file)
@@ -4,11 +4,15 @@
 
 #include "mgraph.h"
 
-/* filled in by characterise_input: */
-static int shift, oldxbits, oldybits, oldx, oldy, oldsz, inc;
+#define OXBITS (XBITS-1)
+#define OX (1<<OXBITS)
+#define OY (Y/2 + 1)
+#define ON (OX*OY)
+#define INC 2
 
 /* filled in by read_input: */
 static struct Vertices all;
+static int computed_count[N];
 
 static void characterise_input(void) {
   struct stat stab;
@@ -18,41 +22,30 @@ static void characterise_input(void) {
 
   if (!stab.st_size || stab.st_size % (sizeof(double)*D3) ||
       stab.st_size > INT_MAX)
-    fail("input file is not reasonable whole number of doubles\n");
-
-  oldsz= stab.st_size / (sizeof(double)*D3);
-  for (shift=1;
-       shift < XBITS+1 && shift < YBITS+1;
-       shift++) {
-    oldxbits= XBITS-1;
-    oldybits= YBITS-1;
-    oldx=  1<<oldxbits;
-    oldy= (1<<oldybits)-1;
-    fprintf(stderr,"sizeof(double)=%d XYBITS=%d,%d, XY=%d*%d=%d"
-           " oldsz=%d shift=%d oldxybits=%d,%d oldxy=%d*%d=%d\n",
-           (int)sizeof(double), XBITS,YBITS, X,Y,N,
-           oldsz,shift,oldxbits,oldybits,oldx,oldy,oldx*oldy);
-    if (oldx*oldy == oldsz) goto found;
-  }
-  fail("input file size cannot be interpolated to target file size\n");
+    fail("input file is not reasonable whole number of vertices\n");
+
+  fprintf(stderr,"XBITS=%d XY=%d*%d=%d\n", XBITS, X,Y,N);
+  fprintf(stderr,"OXBITS=%d OXY=%d*%d=%d\n", OXBITS, OX,OY,ON);
+  fprintf(stderr,"sizeof(double)=%d st_size=%lu\n",
+         (int)sizeof(double), (unsigned long)stab.st_size);
 
- found:
-  inc= 1<<shift;
-  fprintf(stderr,"inc=%d\n",inc);
+  if (stab.st_size != (ON*sizeof(double)*D3))
+    fail("input file wrong size\n");
 }
 
 static void read_input(void) {
   int x,y, ox,oy, v,ov;
   
-  for (oy=y=0; y<Y; oy++, y+=inc) {
+  for (oy=y=0; y<Y; oy++, y+=INC) {
     fprintf(stderr, "y=%2d>%2d", oy,y);
-    for (ox=x=0; x<X; ox++, x+=inc) {
+    for (ox=x=0; x<X; ox++, x+=INC) {
       errno= 0;
-      ov= (oy << oldxbits) | ox;
+      ov= (oy << OXBITS) | ox;
       v= (y << YSHIFT) | x;
       fprintf(stderr, " 0%02o->0x%02x", ov, v);
       if (fread(all.a[v], sizeof(double), D3, stdin) != D3)
        diee("\nread input");
+      computed_count[v]= -1;
     }
     fputc('\n',stderr);
   }
@@ -74,64 +67,116 @@ static void read_input(void) {
    *
    */
 
-#define NEXTV (v= EDGE_END2(v,direction))
-
+typedef struct {
+  int v, e;
+} Traverse;
+
+static void traverse_next(Traverse *t) {
+  int v2;
+  v2= EDGE_END2(t->v, t->e);
+  if (v2>=0) {
+    int e2= edge_reverse(t->v, t->e);
+    assert(EDGE_END2(v2,e2) == t->v);
+    t->e= (e2+3) % V6;
+  }
+  t->v= v2;
+}
+     
 static void interpolate_line(int startvertex,
                             int direction /* edge number */,
-                            int norig) {
-  double xa[norig+1], ya[D3][norig+1], n;
+                            int nmax) {
+  double xa[nmax], ya[D3][nmax];
   const gsl_interp_type *it;
   gsl_interp *interp;
   gsl_interp_accel *accel;
-  int i, v, k;
+  int i, k, nold, nnew;
+  Traverse traverse;
+
+#define STARTV (traverse.v=startvertex, traverse.e=direction)
+#define NEXTV (traverse_next(&traverse), fprintf(stderr,"-%02x",traverse.v))
+
+  fprintf(stderr,"interpolate_line 0x%02x,%d nmax=%d ",
+         startvertex,direction,nmax);
 
-  for (i=0, v=startvertex;
+  for (i=0, STARTV;
        ;
-       i++, NEXTV, assert(v>=0), NEXTV) {
-    if (v<0)
-      break;
-    assert(i <= norig);
+       ) {
+    assert(traverse.v>=0);
+    assert(i < nmax);
     xa[i]= i;
-    K ya[k][i]= all.a[v][k];
-    if (v==startvertex)
-      break;
+    K ya[k][i]= all.a[traverse.v][k];
+    fputc('*',stderr);
+    if (i && traverse.v==startvertex) { fputc('@',stderr); break; }
+    i++;
+    NEXTV;
+    if (traverse.v<0) break;
+    NEXTV;
   }
-  n= i;
-  it= v>=0 ? gsl_interp_cspline_periodic : gsl_interp_cspline;
+  if (traverse.v>=0) {
+    it= gsl_interp_cspline_periodic;
+    nold= i+1;
+    nnew= i;
+  } else {
+    it= gsl_interp_cspline;
+    nold= i;
+    nnew= i-1;
+  }
+
+  fprintf(stderr,"  n=%d->%d loop=0x%02x", nold,nnew,traverse.v);
 
   K {
-    interp= gsl_interp_alloc(it,n); if (!interp) diee("gsl_interp_alloc");
+    fprintf(stderr,"\n  k=%d ", k);
+
+    interp= gsl_interp_alloc(it,nold); 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) );
+    GA( gsl_interp_init(interp,xa,ya[k],nold) );
 
-    for (i=0, v=startvertex;
-        i<n-1;
+    for (i=0, STARTV;
+        i<nnew;
         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]) );
+      assert(traverse.v>=0); NEXTV; assert(traverse.v>=0);
+      fputc('#',stderr);
+      assert(computed_count[traverse.v] == k);
+      GA( gsl_interp_eval_e(interp,xa,ya[k], i+0.5, accel,
+                           &all.a[traverse.v][k]) );
+      computed_count[traverse.v]++;
     }
     
     gsl_interp_accel_free(accel);
     gsl_interp_free(interp);
   }
+  fprintf(stderr,"    done\n");
 }
        
 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 (y=0; y<(Y+1)/2; y+=INC) {
+    interpolate_line(y<<YSHIFT, 0, OX*2+1);
   }
-  for (x=0; x<X; x+=inc) {
-    interpolate_line(x, 5, oldy);
-    interpolate_line(x, 4, oldy);
+  for (x=0; x<X; x+=INC) {
+    interpolate_line(                x, 5, OY);
+    interpolate_line((Y-1)<<YSHIFT | x, 1, OY);
   }
 }
 
 static void write_output(void) {
+  int x, y, v, c, bad=0;
+  for (y=0; y<Y; y++) {
+    fprintf(stderr,"checking y=%d ", y);
+    for (x=0; x<X; x++) {
+      v= y<<YSHIFT | x;
+      fprintf(stderr," %02x",v);
+      c= computed_count[v];
+      if (c==D3) fputc('#',stderr);
+      else if (c==-1) fputc('*',stderr);
+      else { fputc('!',stderr); bad++; }
+    }
+    fputc('\n',stderr);
+  }
+  assert(!bad);
   if (fwrite(all.a,sizeof(all.a),1,stdout) != 1) diee("fwrite");
+  fprintf(stderr,"interpolated\n");
 }
 
 int main(int argc, const char **argv) {
index 4b708e5..7f5f3bd 100644 (file)
--- a/mgraph.c
+++ b/mgraph.c
@@ -11,7 +11,7 @@ static const unsigned dx[2][V6]= {{  +1,   0,  -1,  -1,  -1,   0  },
 int edge_end2(unsigned v1, int e) {
   unsigned x, y;
 
-  y= (v1 & YMASK) + dy[e];
+  y= (v1 & ~XMASK) + dy[e];
   if (y >= Y*Y1) return -1;
 
   x= (v1 & XMASK) + dx[(v1 >> YSHIFT) & 1][e];
@@ -30,9 +30,12 @@ static const unsigned reverse[2][V6]= {{ 3, 4, 5, 0, 1, 2 },
 
 int edge_reverse(int v1, int e) {
   unsigned x2;
-  int flip;
+  int flip, eprime;
 
   x2= (v1 & XMASK) + dx[(v1 >> YSHIFT) & 1][e];
   flip= !!(x2 & ~XMASK);
-  return reverse[flip][e];
+  eprime= reverse[flip][e];
+//  printf("%60s %02x -%d->,<-%d- (%02x) [x2=%x flip=%d]\n","",
+//      v1,e,eprime, EDGE_END2(v1,e), x2,flip);
+  return eprime;
 }
index e98fa85..a5eb611 100644 (file)
--- a/mgraph.h
+++ b/mgraph.h
@@ -60,8 +60,8 @@
  *                  /  \
  *                4/   5\
  *
- * vertex number:   0000 | y     | x
- *                        YBITS   XBITS
+ * vertex number:   0000 | y      | x
+ *                        (YBITS)   XBITS
  */
 
 #ifndef MGRAPH_H
 
 #include "common.h"
 
-#ifndef DEFBITS
+#ifndef DEFSZ /* DEFSZ may be (Y/2-1)*10 + XBITS ie Y is 2*<tens>+1 */
 #define XBITS 3
 #define YBITS 3
+#define Y ((1<<YBITS) - 1)
+#define YMASK (Y << YSHIFT)
 #else
-#define XBITS (DEFBITS / 10)
-#define YBITS (DEFBITS % 10)
+#define XBITS  (DEFSZ % 10)
+#define Y      ((DEFSZ / 10)*2+1)
 #endif
 
 #define X (1<<XBITS)
-#define Y ((1<<YBITS) - 1)
 
 #define N (X*Y)
 #define XMASK (X-1)
 #define YSHIFT XBITS
 #define Y1 (1 << YSHIFT)
-#define YMASK (Y << YSHIFT)
 
 #define V6 6
 
 int edge_end2(unsigned v1, int e);
 #define EDGE_END2 edge_end2
 
-/* given v1,e  s.t.  v2==EDGE_END2(v1,e) >= 0,
- * returns  eprime s.t. v1==EDGE_END2(v2,eprime) */
+/* given    v1,e     s.t.  v2==EDGE_END2(v1,e) >= 0,
+ * returns  eprime   s.t.  v1==EDGE_END2(v2,eprime) */
 int edge_reverse(int v1, int e);
 
-#define RIM_VERTEX_P(v) (((v) & YMASK) == 0 || ((v) & YMASK) == (Y-1)*Y1)
+#define RIM_VERTEX_P(v) (((v) & ~XMASK) == 0 || ((v) & ~XMASK) == (Y-1)*Y1)
 
 #define FOR_VEDGE_X(v1,e,v2,init,otherwise)    \
   FOR_VPEDGE((v1),(e))                         \