From 75293637fe7efbce3f2e873cc792da3bff577fef Mon Sep 17 00:00:00 2001 From: Ian Jackson Date: Fri, 18 Jan 2008 23:22:44 +0000 Subject: [PATCH] interpolates but not very well --- .bzrignore | 4 +- Makefile | 10 ++-- interpolate.c | 151 ++++++++++++++++++++++++++++++++------------------ mgraph.c | 9 ++- mgraph.h | 20 +++---- 5 files changed, 121 insertions(+), 73 deletions(-) diff --git a/.bzrignore b/.bzrignore index 93c99e1..156d40a 100644 --- a/.bzrignore +++ b/.bzrignore @@ -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 diff --git a/Makefile b/Makefile index a2409fc..dd4b218 100644 --- 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 diff --git a/interpolate.c b/interpolate.c index d9d3910..1de2bea 100644 --- a/interpolate.c +++ b/interpolate.c @@ -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< 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<%2d", oy,y); - for (ox=x=0; x0x%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=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*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; } diff --git a/mgraph.h b/mgraph.h index e98fa85..a5eb611 100644 --- 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 @@ -69,22 +69,22 @@ #include "common.h" -#ifndef DEFBITS +#ifndef DEFSZ /* DEFSZ may be (Y/2-1)*10 + XBITS ie Y is 2*+1 */ #define XBITS 3 #define YBITS 3 +#define Y ((1<= 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)) \ -- 2.30.2