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++) {
63 if (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]= -1;
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 const gsl_interp_type *it;
178 gsl_interp_accel *accel;
179 int i, k, nold, nnew;
182 #define STARTV (traverse.v=startvertex, traverse.e=direction)
183 #define NEXTV (traverse_next(&traverse), printf("-%02x",traverse.v))
185 printf("interpolate_line 0x%02x,%d nmax=%d ",
186 startvertex,direction,nmax);
191 assert(traverse.v>=0);
194 K ya[k][i]= all.a[traverse.v][k];
196 if (i && traverse.v==startvertex) { putchar('@'); break; }
199 if (traverse.v<0) break;
212 printf(" n=%d->%d loop=0x%02x", nold,nnew,traverse.v);
215 printf("\n k=%d ", k);
217 interp= gsl_interp_alloc(it,nold); if (!interp) diee("gsl_interp_alloc");
218 accel= gsl_interp_accel_alloc(); if(!accel) diee("gsl_interp_accel_alloc");
219 GA( gsl_interp_init(interp,xa,ya[k],nold) );
224 assert(traverse.v>=0); NEXTV; assert(traverse.v>=0);
226 GA( gsl_interp_eval_e(interp,xa,ya[k], i+0.5, accel,
227 &all.a[traverse.v][k]) );
228 note_computed(traverse.v, k);
231 gsl_interp_accel_free(accel);
232 gsl_interp_free(interp);
237 static void interpolate_gsl(const gsl_interp_type *aperiodic,
238 const gsl_interp_type *periodic) {
241 for (y=0; y<(Y+1)/2; y+=inc) {
242 interpolate_gsl_line(y<<YSHIFT | ((y>>1)&1), 0, oX*2+1,
243 aperiodic, periodic);
245 for (x=0; x<X; x+=inc) {
246 interpolate_gsl_line( x, 5, oY, aperiodic, periodic);
247 interpolate_gsl_line((Y-1)<<YSHIFT | x, 1, oY, aperiodic, periodic);
251 static void write_output(const char *fn) {
254 int x, y, v, c, bad=0;
256 for (y=0; y<Y; y++) {
257 printf("checking y=%d ", y);
258 for (x=0; x<X; x++) {
261 c= computed_count[v];
262 if (c==D3) putchar('#');
263 else if (c==-1) putchar('*');
264 else { printf("!%d",c); bad++; }
270 if (asprintf(&fn_new,"%s.new",fn) <= 0) diee("asprintf");
271 f= fopen(fn_new,"wb"); if (!f) diee("open new output file");
272 if (fwrite(all.a,sizeof(all.a),1,f) != 1) diee("fwrite");
273 if (fclose(f)) diee("fclose output");
274 if (rename(fn_new,fn)) diee("install new output");
275 printf("interpolated\n");
278 int main(int argc, const char **argv) {
280 strncmp(argv[1],"-a",2) || strlen(argv[1])!=3 ||
282 strncmp(argv[3],"-o",2)) {
283 fputs("usage: interpolate -aK INPUT.cfm -oOUTPUT.cfm\n",stderr);
287 if (!freopen(argv[2],"rb",stdin)) diee("open input");
288 characterise_input();
291 switch (argv[1][2]) {
292 case '4': interpolate_lin4point(); break;
293 case 'a': interpolate_gsl(gsl_interp_akima,
294 gsl_interp_akima_periodic); break;
295 case 'c': interpolate_gsl(gsl_interp_cspline,
296 gsl_interp_cspline_periodic); break;
298 fail("unknown interpolation algorithm\n");
301 write_output(argv[3]+2);