2 * increases the resolution of a conformation by interpolation
7 /* filled in by characterise_input: */
8 static int oXBITS,oX,oY,oN,inc;
10 /* filled in by read_input: */
11 static struct Vertices all;
12 static int computed_count[N];
14 static void note_computed(int v, int k) {
15 assert(computed_count[v] == k);
19 static void characterise_input(void) {
21 int r, shift, oN_calc;
23 r= fstat(0,&stab); if (r) diee("fstat input to find length");
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");
29 printf("XBITS=%d XY=%d*%d=%d", XBITS, X,Y,N);
31 oN= stab.st_size / (sizeof(double)*D3);
36 printf("\n shift=%d inc=%d ",shift,inc);
38 if (X % inc) continue;
39 oX= X/inc; printf("oX=%d ",oX);
41 if ((Y-1) % inc) continue;
42 oY= (Y-1)/inc + 1; printf("oY=%d ",oY);
45 printf("oN=%d",oN_calc);
46 if (oN_calc == oN) goto found;
48 fail("\ninput size cannot be reconciled\n");
52 printf("oXBITS=%d\n", oXBITS);
55 static void read_input(void) {
58 for (oy=0; oy<oY; oy++) {
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; }
65 ov= (oy << oXBITS) | ox;
67 printf(" 0%02o->0x%02x", ov, v);
68 if (fread(all.a[v], sizeof(double), D3, stdin) != D3)
70 computed_count[v]= D3*2;
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.
82 * This directly gives the required additional vertices:
84 * ,O. O original vertex
86 * ,* - - *. dashed lines are new edges,
87 * / ` ' \ not directly represented
88 * O----`*'----O or manipulated by this program
96 static void traverse_next(Traverse *t) {
102 v2= EDGE_END2(t->v, t->e);
104 int e2= edge_reverse(t->v, t->e);
105 assert(EDGE_END2(v2,e2) == t->v);
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];
117 if (inc != 2) fail("cannot do >2x lin4point interpolation\n");
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;
127 traverse_next(&trav);
131 traverse_next(&trav);
135 traverse_next(&trav);
136 traverse_next(&trav);
139 trav.v= vq; trav.e= EDGE_OPPOSITE(eqr);
140 traverse_next(&trav);
141 traverse_next(&trav);
144 printf(" 0x%02x-##-%02x-!%02x!-%02x-##-%02x",
145 vp&0xff,vq,vm,vr,vs&0xff);
148 K pqtarg[k]= all.a[vq][k]*2 - all.a[vp][k];
150 K pqtarg[k]= all.a[vr][k];
153 K srtarg[k]= all.a[vr][k]*2 - all.a[vs][k];
155 K srtarg[k]= all.a[vq][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])
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;
179 gsl_interp_accel *accel;
180 int i, k, nold, nnew, x, skipped, realstart;
183 printf("interpolate_line 0x%02x,%d nmax=%d ",
184 startvertex,direction,nmax);
186 traverse.v=startvertex, traverse.e=direction;
188 while (!computed_count[traverse.v]) {
189 printf("-%02x~",traverse.v);
190 traverse_next(&traverse);
192 assert(skipped < nmax);
194 realstart= traverse.v;
199 x++, traverse_next(&traverse)) {
200 printf("-%02x",traverse.v);
207 if (computed_count[traverse.v]) {
210 K ya[k][nold]= all.a[traverse.v][k];
216 vnew[nnew]= traverse.v;
220 if (x && traverse.v==realstart) {
227 printf(" n=%d->%d", nold,nnew);
228 if (!nnew) { printf(" nothing\n"); return; }
231 printf(" x=%g..%g->%d..%d", xa[0],xa[nold-1], xnew[0],xnew[nnew-1]);
234 printf("\n k=%d ", k);
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) );
240 for (i=0; i<nnew; i++) {
243 printf(" %02x(%d)",v,x);
244 GA( gsl_interp_eval_e(interp,xa,ya[k], x, accel, &all.a[v][k]) );
248 gsl_interp_accel_free(accel);
249 gsl_interp_free(interp);
254 static void interpolate_gsl(const gsl_interp_type *aperiodic,
255 const gsl_interp_type *periodic) {
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);
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);
268 static void write_output(const char *fn) {
271 int x, y, v, c, bad=0;
273 for (y=0; y<Y; y++) {
274 printf("checking y=%d ", y);
275 for (x=0; x<X; x++) {
278 c= computed_count[v];
279 if (c==D3) putchar('#');
280 else if (c==D3*2) putchar('*');
281 else { printf("!%d",c); bad++; }
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");
295 int main(int argc, const char **argv) {
297 strncmp(argv[1],"-a",2) || strlen(argv[1])!=3 ||
299 strncmp(argv[3],"-o",2)) {
300 fputs("usage: interpolate -aK INPUT.cfm -oOUTPUT.cfm\n",stderr);
304 if (!freopen(argv[2],"rb",stdin)) diee("open input");
307 characterise_input();
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;
317 fail("unknown interpolation algorithm\n");
320 write_output(argv[3]+2);