#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");
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');
}
}
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);
}
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_cspline_periodic;
- nold= i+1;
- nnew= i;
- } else {
- it= gsl_interp_cspline;
- 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");
+ mgraph_prepare();
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;
}
~ian / moebius2.git / blobdiff