interpolate.c

   1 /*
   2  * increases the resolution of a conformation by interpolation
   3  */
   4
   5 #include "mgraph.h"
   6
   7 /* filled in by characterise_input: */
   8 static int oXBITS,oX,oY,oN,inc;
   9
  10 /* filled in by read_input: */
  11 static struct Vertices all;
  12 static int computed_count[N];
  13
  14 static void note_computed(int v, int k) {
  15   assert(computed_count[v] == k);
  16   computed_count[v]++;
  17 }
  18
  19 static void characterise_input(void) {
  20   struct stat stab;
  21   int r, shift, oN_calc;
  22
  23   r= fstat(0,&stab);  if (r) diee("fstat input to find length");
  24
  25   if (!stab.st_size || stab.st_size % (sizeof(double)*D3) ||
  26       stab.st_size > INT_MAX)
  27     fail("input file is not reasonable whole number of vertices\n");
  28
  29   printf("XBITS=%d XY=%d*%d=%d", XBITS, X,Y,N);
  30
  31   oN= stab.st_size / (sizeof(double)*D3);
  32
  33   for (shift=1, inc=2;
  34        shift<XBITS-1;
  35        shift++, inc<<=1) {
  36     printf("\n shift=%d inc=%d ",shift,inc);
  37
  38     if (X % inc) continue;
  39     oX= X/inc; printf("oX=%d ",oX);
  40
  41     if ((Y-1) % inc) continue;
  42     oY= (Y-1)/inc + 1; printf("oY=%d ",oY);
  43
  44     oN_calc= oX*oY;
  45     printf("oN=%d",oN_calc);
  46     if (oN_calc == oN) goto found;
  47   }
  48   fail("\ninput size cannot be reconciled\n");
  49
  50  found:
  51   oXBITS= XBITS-shift;
  52   printf("oXBITS=%d\n", oXBITS);
  53 }
  54
  55 static void read_input(void) {
  56   int x,y, ox,oy, v,ov;
  57
  58   for (oy=0; oy<oY; oy++) {
  59     y= oy*inc;
  60     printf("y=%2d>%2d", oy,y);
  61     for (ox=0; ox<oX; ox++) {
  62       x= (ox - (oy>>1))*inc + (y>>1) + X*2;
  63       while (x>=X) { y= (Y-1)-y; x-=X; }
  64       errno= 0;
  65       ov= (oy << oXBITS) | ox;
  66       v= (y << YSHIFT) | x;
  67       printf(" 0%02o->0x%02x", ov, v);
  68       if (fread(all.a[v], sizeof(double), D3, stdin) != D3)
  69         diee("\nread input");
  70       computed_count[v]= D3*2;
  71     }
  72     putchar('\n');
  73   }
  74 }
  75
  76   /* We use GSL's interpolation functions.  Each xa array is simple
  77    * integers, and we do the interpolation three times for each line
  78    * of topologically colinear vertices (once, separately, for each of
  79    * the three spatial coordinates in the ya array).  The `horizontal'
  80    * lines are done with periodic boundary conditions, obviously.
  81    *
  82    * This directly gives the required additional vertices:
  83    *
  84    *         ,O.                   O  original vertex
  85    *        /   \                  *  new vertex
  86    *      ,* - - *.                dashed lines are new edges,
  87    *     /  `   '  \                not directly represented
  88    *    O----`*'----O               or manipulated by this program
  89    *
  90    */
  91
  92 typedef struct {
  93   int v, e;
  94 } Traverse;
  95
  96 static void traverse_next(Traverse *t) {
  97   int v2;
  98
  99   if (t->v<0)
 100     return;
 101
 102   v2= EDGE_END2(t->v, t->e);
 103   if (v2>=0) {
 104     int e2= edge_reverse(t->v, t->e);
 105     assert(EDGE_END2(v2,e2) == t->v);
 106     t->e= (e2+3) % V6;
 107   }
 108   t->v= v2;
 109 }
 110
 111 static void interpolate_lin4point(void) {
 112   /* four points P Q R S, although P and S may be missing
 113    * interpolate in QR finding M. */
 114   int xq,yq,eqr, vp,vq,vm,vr,vs, k;
 115   double pqtarg[D3], srtarg[D3];
 116
 117   if (inc != 2) fail("cannot do >2x lin4point interpolation\n");
 118
 119   for (eqr=1; eqr!=4; eqr+=5, eqr%=V6) { /* each old edge exactly once */
 120     printf("eqr=%d\n",eqr);
 121     for (yq=0; yq<Y; yq+=inc) {
 122       printf(" yq=%2d ",yq);
 123       for (xq=((yq>>1)&1); xq<X; xq+=inc) {
 124         vq= yq << YSHIFT | xq;
 125         Traverse trav; trav.v=vq; trav.e=eqr;
 126
 127         traverse_next(&trav);
 128         vm= trav.v;
 129         if (vm<0) continue;
 130
 131         traverse_next(&trav);
 132         vr= trav.v;
 133         assert(vr>=0);
 134
 135         traverse_next(&trav);
 136         traverse_next(&trav);
 137         vs= trav.v;
 138
 139         trav.v= vq; trav.e= EDGE_OPPOSITE(eqr);
 140         traverse_next(&trav);
 141         traverse_next(&trav);
 142         vp= trav.v;
 143
 144         printf(" 0x%02x-##-%02x-!%02x!-%02x-##-%02x",
 145                 vp&0xff,vq,vm,vr,vs&0xff);
 146
 147         if (vp>=0)
 148           K pqtarg[k]= all.a[vq][k]*2 - all.a[vp][k];
 149         else
 150           K pqtarg[k]= all.a[vr][k];
 151
 152         if (vs>=0)
 153           K srtarg[k]= all.a[vr][k]*2 - all.a[vs][k];
 154         else
 155           K srtarg[k]= all.a[vq][k];
 156
 157         K {
 158           const double alpha= 6;
 159           all.a[vm][k]= 0.5 * (all.a[vq][k] + all.a[vr][k]);
 160           pqtarg[k]= 0.5 * (pqtarg[k] + all.a[vm][k]);
 161           srtarg[k]= 0.5 * (srtarg[k] + all.a[vm][k]);
 162           all.a[vm][k]= (pqtarg[k] + srtarg[k] + alpha * all.a[vm][k])
 163             / (2 + alpha);
 164           note_computed(vm,k);
 165         }
 166       }
 167       putchar('\n');
 168     }
 169   }
 170 }
 171
 172 static void interpolate_gsl_line
 173   (int startvertex, int direction /* edge number */, int nmax,
 174    const gsl_interp_type *aperiodic, const gsl_interp_type *periodic) {
 175   double xa[nmax], ya[D3][nmax];
 176   int vnew[nmax], xnew[nmax];
 177   const gsl_interp_type *it;
 178   gsl_interp *interp;
 179   gsl_interp_accel *accel;
 180   int i, k, nold, nnew, x, skipped, realstart;
 181   Traverse traverse;
 182
 183   printf("interpolate_line 0x%02x,%d nmax=%d ",
 184           startvertex,direction,nmax);
 185
 186   traverse.v=startvertex, traverse.e=direction;
 187   skipped= 0;
 188   while (!computed_count[traverse.v]) {
 189     printf("-%02x~",traverse.v);
 190     traverse_next(&traverse);
 191     skipped++;
 192     assert(skipped < nmax);
 193   }
 194   realstart= traverse.v;
 195
 196   nold=0, nnew=0;
 197   for (x=0;
 198        ;
 199        x++, traverse_next(&traverse)) {
 200     printf("-%02x",traverse.v);
 201     if (traverse.v<0) {
 202       putchar('$');
 203       it= aperiodic;
 204       assert(!skipped);
 205       break;
 206     }
 207     if (computed_count[traverse.v]) {
 208       assert(nold < nmax);
 209       xa[nold]= x;
 210       K ya[k][nold]= all.a[traverse.v][k];
 211       nold++;
 212       putchar('*');
 213     } else {
 214       assert(nnew < nmax);
 215       xnew[nnew]= x;
 216       vnew[nnew]= traverse.v;
 217       nnew++;
 218       putchar('#');
 219     }
 220     if (x && traverse.v==realstart) {
 221       putchar('@');
 222       it= periodic;
 223       break;
 224     }
 225   }
 226
 227   printf("  n=%d->%d", nold,nnew);
 228   if (!nnew) { printf(" nothing\n"); return; }
 229
 230   assert(nold);
 231   printf(" x=%g..%g->%d..%d", xa[0],xa[nold-1], xnew[0],xnew[nnew-1]);
 232
 233   K {
 234     printf("\n  k=%d ", k);
 235
 236     interp= gsl_interp_alloc(it,nold); if (!interp) diee("gsl_interp_alloc");
 237     accel= gsl_interp_accel_alloc(); if(!accel) diee("gsl_interp_accel_alloc");
 238     GA( gsl_interp_init(interp,xa,ya[k],nold) );
 239
 240     for (i=0; i<nnew; i++) {
 241       int v= vnew[i];
 242       x= xnew[i];
 243       printf(" %02x(%d)",v,x);
 244       GA( gsl_interp_eval_e(interp,xa,ya[k], x, accel, &all.a[v][k]) );
 245       note_computed(v,k);
 246     }
 247
 248     gsl_interp_accel_free(accel);
 249     gsl_interp_free(interp);
 250   }
 251   printf("    done\n");
 252 }
 253
 254 static void interpolate_gsl(const gsl_interp_type *aperiodic,
 255                             const gsl_interp_type *periodic) {
 256   int x,y;
 257
 258   for (x=0; x<X; x+=inc) {
 259     interpolate_gsl_line(                x, 5, Y, aperiodic, periodic);
 260     interpolate_gsl_line((Y-1)<<YSHIFT | x, 1, Y, aperiodic, periodic);
 261   }
 262   for (y=0; y<(Y+1)/2; y++) {
 263     interpolate_gsl_line(y<<YSHIFT | ((y>>1)&1), 0, X*2+1,
 264                          aperiodic, periodic);
 265   }
 266 }
 267
 268 static void write_output(const char *fn) {
 269   FILE *f;
 270   char *fn_new;
 271   int x, y, v, c, bad=0;
 272
 273   for (y=0; y<Y; y++) {
 274     printf("checking y=%d ", y);
 275     for (x=0; x<X; x++) {
 276       v= y<<YSHIFT | x;
 277       printf(" %02x",v);
 278       c= computed_count[v];
 279       if (c==D3) putchar('#');
 280       else if (c==D3*2) putchar('*');
 281       else { printf("!%d",c); bad++; }
 282     }
 283     putchar('\n');
 284   }
 285   assert(!bad);
 286
 287   if (asprintf(&fn_new,"%s.new",fn) <= 0) diee("asprintf");
 288   f= fopen(fn_new,"wb");  if (!f) diee("open new output file");
 289   if (fwrite(all.a,sizeof(all.a),1,f) != 1) diee("fwrite");
 290   if (fclose(f)) diee("fclose output");
 291   if (rename(fn_new,fn)) diee("install new output");
 292   printf("interpolated\n");
 293 }
 294
 295 int main(int argc, const char **argv) {
 296   if (argc!=4 ||
 297       strncmp(argv[1],"-a",2) || strlen(argv[1])!=3 ||
 298       argv[2][0]=='-' ||
 299       strncmp(argv[3],"-o",2)) {
 300     fputs("usage: interpolate -aK INPUT.cfm -oOUTPUT.cfm\n",stderr);
 301     exit(8);
 302   }
 303
 304   if (!freopen(argv[2],"rb",stdin)) diee("open input");
 305
 306   mgraph_prepare();
 307   characterise_input();
 308   read_input();
 309
 310   switch (argv[1][2]) {
 311   case '4':  interpolate_lin4point();                       break;
 312   case 'a':  interpolate_gsl(gsl_interp_akima,
 313                              gsl_interp_akima_periodic);    break;
 314   case 'c':  interpolate_gsl(gsl_interp_cspline,
 315                              gsl_interp_cspline_periodic);  break;
 316   default:
 317     fail("unknown interpolation algorithm\n");
 318   }
 319
 320   write_output(argv[3]+2);
 321   return 0;
 322 }