X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?p=moebius2.git;a=blobdiff_plain;f=interpolate.c;h=a2d1cafd48de65ad50b943f92afd6495a11de0c7;hp=f5d7c1bb9c727517e9bb3784f3cc5f64ac41d9b0;hb=ad8f9959ed7e147b88a57a7634c78faa69cefc6a;hpb=df4fafe44b168dc0e79489e558fc00b3a2b7da3f

diff --git a/interpolate.c b/interpolate.c
index f5d7c1b..a2d1caf 100644
--- a/interpolate.c
+++ b/interpolate.c
@@ -4,19 +4,21 @@
 
 #include "mgraph.h"
 
-#define OXBITS (XBITS-1)
-#define OX (1<<OXBITS)
-#define OY (Y/2 + 1)
-#define ON (OX*OY)
-#define INC 2
+/* filled in by characterise_input: */
+static int oXBITS,oX,oY,oN,inc;
 
 /* filled in by read_input: */
 static struct Vertices all;
 static int computed_count[N];
 
+static void note_computed(int v, int k) {
+  assert(computed_count[v] == k);
+  computed_count[v]++;
+}
+
 static void characterise_input(void) {
   struct stat stab;
-  int r;
+  int r, shift, oN_calc;
   
   r= fstat(0,&stab);  if (r) diee("fstat input to find length");
 
@@ -24,30 +26,50 @@ static void characterise_input(void) {
       stab.st_size > INT_MAX)
     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);
+  printf("XBITS=%d XY=%d*%d=%d", XBITS, X,Y,N);
+
+  oN= stab.st_size / (sizeof(double)*D3);
+
+  for (shift=1, inc=2;
+       shift<XBITS-1;
+       shift++, inc<<=1) {
+    printf("\n shift=%d inc=%d ",shift,inc);
+
+    if (X % inc) continue;
+    oX= X/inc; printf("oX=%d ",oX);
 
-  if (stab.st_size != (ON*sizeof(double)*D3))
-    fail("input file wrong size\n");
+    if ((Y-1) % inc) continue;
+    oY= (Y-1)/inc + 1; printf("oY=%d ",oY);
+
+    oN_calc= oX*oY;
+    printf("oN=%d",oN_calc);
+    if (oN_calc == oN) goto found;
+  }
+  fail("\ninput size cannot be reconciled\n");
+
+ found:
+  oXBITS= XBITS-shift;
+  printf("oXBITS=%d\n", oXBITS);
 }
 
 static void read_input(void) {
   int x,y, ox,oy, v,ov;
   
-  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 (oy=0; oy<oY; oy++) {
+    y= oy*inc;
+    printf("y=%2d>%2d", oy,y);
+    for (ox=0; ox<oX; ox++) {
+      x= (ox - (oy>>1))*inc + (y>>1) + X*2;
+      while (x>=X) { y= (Y-1)-y; x-=X; }
       errno= 0;
-      ov= (oy << OXBITS) | ox;
+      ov= (oy << oXBITS) | ox;
       v= (y << YSHIFT) | x;
-      fprintf(stderr, " 0%02o->0x%02x", ov, v);
+      printf(" 0%02o->0x%02x", ov, v);
       if (fread(all.a[v], sizeof(double), D3, stdin) != D3)
 	diee("\nread input");
-      computed_count[v]= -1;
+      computed_count[v]= D3*2;
     }
-    fputc('\n',stderr);
+    putchar('\n');
   }
 }
 
@@ -73,6 +95,10 @@ typedef struct {
 
 static void traverse_next(Traverse *t) {
   int v2;
+
+  if (t->v<0)
+    return;
+  
   v2= EDGE_END2(t->v, t->e);
   if (v2>=0) {
     int e2= edge_reverse(t->v, t->e);
@@ -81,112 +107,214 @@ static void traverse_next(Traverse *t) {
   }
   t->v= v2;
 }
-     
-static void interpolate_line(int startvertex,
-			     int direction /* edge number */,
-			     int nmax) {
+
+static void interpolate_lin4point(void) {
+  /* four points P Q R S, although P and S may be missing
+   * interpolate in QR finding M. */
+  int xq,yq,eqr, vp,vq,vm,vr,vs, k;
+  double pqtarg[D3], srtarg[D3];
+
+  if (inc != 2) fail("cannot do >2x lin4point interpolation\n");
+  
+  for (eqr=1; eqr!=4; eqr+=5, eqr%=V6) { /* each old edge exactly once */
+    printf("eqr=%d\n",eqr);
+    for (yq=0; yq<Y; yq+=inc) {
+      printf(" yq=%2d ",yq);
+      for (xq=((yq>>1)&1); xq<X; xq+=inc) {
+	vq= yq << YSHIFT | xq;
+	Traverse trav; trav.v=vq; trav.e=eqr;
+
+	traverse_next(&trav);
+	vm= trav.v;
+	if (vm<0) continue;
+
+	traverse_next(&trav);
+	vr= trav.v;
+	assert(vr>=0);
+
+	traverse_next(&trav);
+	traverse_next(&trav);
+	vs= trav.v;
+
+	trav.v= vq; trav.e= EDGE_OPPOSITE(eqr);
+	traverse_next(&trav);
+	traverse_next(&trav);
+	vp= trav.v;
+	
+	printf(" 0x%02x-##-%02x-!%02x!-%02x-##-%02x",
+		vp&0xff,vq,vm,vr,vs&0xff);
+
+	if (vp>=0)
+	  K pqtarg[k]= all.a[vq][k]*2 - all.a[vp][k];
+	else
+	  K pqtarg[k]= all.a[vr][k];
+
+	if (vs>=0)
+	  K srtarg[k]= all.a[vr][k]*2 - all.a[vs][k];
+	else
+	  K srtarg[k]= all.a[vq][k];
+
+	K {
+	  const double alpha= 6;
+	  all.a[vm][k]= 0.5 * (all.a[vq][k] + all.a[vr][k]);
+	  pqtarg[k]= 0.5 * (pqtarg[k] + all.a[vm][k]);
+	  srtarg[k]= 0.5 * (srtarg[k] + all.a[vm][k]);
+	  all.a[vm][k]= (pqtarg[k] + srtarg[k] + alpha * all.a[vm][k])
+	    / (2 + alpha);
+	  note_computed(vm,k);
+	}
+      }
+      putchar('\n');
+    }
+  }
+}
+
+static void interpolate_gsl_line
+  (int startvertex, int direction /* edge number */, int nmax,
+   const gsl_interp_type *aperiodic, const gsl_interp_type *periodic) {
   double xa[nmax], ya[D3][nmax];
+  int vnew[nmax], xnew[nmax];
   const gsl_interp_type *it;
   gsl_interp *interp;
   gsl_interp_accel *accel;
-  int i, k, nold, nnew;
+  int i, k, nold, nnew, x, skipped, realstart;
   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 ",
+  printf("interpolate_line 0x%02x,%d nmax=%d ",
 	  startvertex,direction,nmax);
 
-  for (i=0, STARTV;
-       ;
-       ) {
-    assert(traverse.v>=0);
-    assert(i < nmax);
-    xa[i]= i;
-    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;
+  traverse.v=startvertex, traverse.e=direction;
+  skipped= 0;
+  while (!computed_count[traverse.v]) {
+    printf("-%02x~",traverse.v);
+    traverse_next(&traverse);
+    skipped++;
+    assert(skipped < nmax);
   }
-  if (traverse.v>=0) {
-    it= gsl_interp_akima_periodic;
-    nold= i+1;
-    nnew= i;
-  } else {
-    it= gsl_interp_akima;
-    nold= i;
-    nnew= i-1;
+  realstart= traverse.v;
+
+  nold=0, nnew=0;
+  for (x=0;
+       ;
+       x++, traverse_next(&traverse)) {
+    printf("-%02x",traverse.v);
+    if (traverse.v<0) {
+      putchar('$');
+      it= aperiodic;
+      assert(!skipped);
+      break;
+    }
+    if (computed_count[traverse.v]) {
+      assert(nold < nmax);
+      xa[nold]= x;
+      K ya[k][nold]= all.a[traverse.v][k];
+      nold++;
+      putchar('*');
+    } else {
+      assert(nnew < nmax);
+      xnew[nnew]= x;
+      vnew[nnew]= traverse.v;
+      nnew++;
+      putchar('#');
+    }
+    if (x && traverse.v==realstart) {
+      putchar('@');
+      it= periodic;
+      break;
+    }
   }
 
-  fprintf(stderr,"  n=%d->%d loop=0x%02x", nold,nnew,traverse.v);
+  printf("  n=%d->%d", nold,nnew);
+  if (!nnew) { printf(" nothing\n"); return; }
+
+  assert(nold);
+  printf(" x=%g..%g->%d..%d", xa[0],xa[nold-1], xnew[0],xnew[nnew-1]);
 
   K {
-    fprintf(stderr,"\n  k=%d ", k);
+    printf("\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],nold) );
 
-    for (i=0, STARTV;
-	 i<nnew;
-	 i++, NEXTV) {
-      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]++;
+    for (i=0; i<nnew; i++) {
+      int v= vnew[i];
+      x= xnew[i];
+      printf(" %02x(%d)",v,x);
+      GA( gsl_interp_eval_e(interp,xa,ya[k], x, accel, &all.a[v][k]) );
+      note_computed(v,k);
     }
     
     gsl_interp_accel_free(accel);
     gsl_interp_free(interp);
   }
-  fprintf(stderr,"    done\n");
+  printf("    done\n");
 }
 	
-static void interpolate(void) {
+static void interpolate_gsl(const gsl_interp_type *aperiodic,
+			    const gsl_interp_type *periodic) {
   int x,y;
   
-  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_gsl_line(                x, 5, Y, aperiodic, periodic);
+    interpolate_gsl_line((Y-1)<<YSHIFT | x, 1, Y, aperiodic, periodic);
   }
-  for (x=0; x<X; x+=INC) {
-    interpolate_line(                x, 5, OY);
-    interpolate_line((Y-1)<<YSHIFT | x, 1, OY);
+  for (y=0; y<(Y+1)/2; y++) {
+    interpolate_gsl_line(y<<YSHIFT | ((y>>1)&1), 0, X*2+1,
+			 aperiodic, periodic);
   }
 }
 
-static void write_output(void) {
+static void write_output(const char *fn) {
+  FILE *f;
+  char *fn_new;
   int x, y, v, c, bad=0;
+  
   for (y=0; y<Y; y++) {
-    fprintf(stderr,"checking y=%d ", y);
+    printf("checking y=%d ", y);
     for (x=0; x<X; x++) {
       v= y<<YSHIFT | x;
-      fprintf(stderr," %02x",v);
+      printf(" %02x",v);
       c= computed_count[v];
-      if (c==D3) fputc('#',stderr);
-      else if (c==-1) fputc('*',stderr);
-      else { fputc('!',stderr); bad++; }
+      if (c==D3) putchar('#');
+      else if (c==D3*2) putchar('*');
+      else { printf("!%d",c); bad++; }
     }
-    fputc('\n',stderr);
+    putchar('\n');
   }
   assert(!bad);
-  if (fwrite(all.a,sizeof(all.a),1,stdout) != 1) diee("fwrite");
-  fprintf(stderr,"interpolated\n");
+
+  if (asprintf(&fn_new,"%s.new",fn) <= 0) diee("asprintf");
+  f= fopen(fn_new,"wb");  if (!f) diee("open new output file");
+  if (fwrite(all.a,sizeof(all.a),1,f) != 1) diee("fwrite");
+  if (fclose(f)) diee("fclose output");
+  if (rename(fn_new,fn)) diee("install new output");
+  printf("interpolated\n");
 }
 
 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); }
+  if (argc!=4 ||
+      strncmp(argv[1],"-a",2) || strlen(argv[1])!=3 ||
+      argv[2][0]=='-' ||
+      strncmp(argv[3],"-o",2)) {
+    fputs("usage: interpolate -aK INPUT.cfm -oOUTPUT.cfm\n",stderr);
+    exit(8);
+  }
 
+  if (!freopen(argv[2],"rb",stdin)) diee("open input");
   characterise_input();
   read_input();
-  interpolate();
-  write_output();
-  if (fclose(stdout)) diee("fclose stdout");
+
+  switch (argv[1][2]) {
+  case '4':  interpolate_lin4point();                       break;
+  case 'a':  interpolate_gsl(gsl_interp_akima,
+			     gsl_interp_akima_periodic);    break;
+  case 'c':  interpolate_gsl(gsl_interp_cspline,
+			     gsl_interp_cspline_periodic);  break;
+  default:
+    fail("unknown interpolation algorithm\n");
+  }
+  
+  write_output(argv[3]+2);
   return 0;
 }