Back to index

wims  3.65+svn20090927
shortpath.c
Go to the documentation of this file.
00001 /*    Copyright (C) 1998-2003 XIAO, Gang of Universite de Nice - Sophia Antipolis
00002  *
00003  *  This program is free software; you can redistribute it and/or modify
00004  *  it under the terms of the GNU General Public License as published by
00005  *  the Free Software Foundation; either version 2 of the License, or
00006  *  (at your option) any later version.
00007  *
00008  *  This program is distributed in the hope that it will be useful,
00009  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
00010  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00011  *  GNU General Public License for more details.
00012  *
00013  *  You should have received a copy of the GNU General Public License
00014  *  along with this program; if not, write to the Free Software
00015  *  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00016  */
00017 
00018  /* Finds the shortest paths linking given points. */
00019 
00020        /* fast algorithm is not ready. */
00021 
00022 #define MAX_POINTS 32       /* cannot exceed number of bits in long */
00023 #define MAX_DIM 10   /* dimension limit */
00024 #define MAX_SHORT 20 /* limit of registered shortest paths */
00025 
00026 #include "../config.h"
00027 #include <stdio.h>
00028 #include <stdlib.h>
00029 #include <string.h>
00030 #include <math.h>
00031 #include <ctype.h>
00032 
00033 int  pointcnt,dimension,first,last,count,shortcnt;
00034        /* style = 0: loop to the start
00035         *        1: arbitrary open path
00036         *        2: open path with fixed start
00037         *        3: open path with fixed end
00038         *        4: open path with fixed start and end */
00039 int style=0;
00040 int fast=0;   /* fast algo if 1, no rigor */
00041 int optimize=0;
00042 char *input;
00043 double dist[MAX_POINTS][MAX_POINTS];
00044 double ds[MAX_POINTS];
00045 double coord[MAX_POINTS][MAX_DIM];
00046 double shortest;
00047 double tolerance; /* to compensate floating rounding error in comparison */
00048 double distmin;      /* minimum of distance between two points */
00049 double firstlast; /* distance between first and last */
00050 int shpath[MAX_SHORT][MAX_POINTS];
00051 int chain[MAX_POINTS];
00052 
00053        /* reverse a subchain */
00054 void reverse(int i, int j)
00055 {
00056     int k,t;
00057     for(k=0;k<(j-i+1)/2;k++) {
00058        t=chain[i+k]; chain[i+k]=chain[j-k];
00059        chain[j-k]=t;
00060     }
00061 }
00062 
00063        /* move point i to position j */
00064 void reinsert(int i, int j)
00065 {
00066     int t;
00067     t=chain[i];
00068     if(j>i) {
00069        j--;
00070        memmove(chain+i,chain+i+1,(j-i)*sizeof(chain[0]));
00071        chain[j]=t;
00072     }
00073     else {
00074        memmove(chain+j+1,chain+j,(i-j)*sizeof(chain[0]));
00075        chain[j]=t;
00076     }
00077 }
00078 
00079 int _fcalc(void)
00080 {
00081     int modif;
00082     int i,j;
00083     
00084     modif=0;
00085        /* four-configuration optimize */
00086     for(i=1;i<pointcnt-2;i++) {
00087        for(j=i+1;j<pointcnt-1;j++) {
00088            if(dist[chain[i-1]][chain[i]]+dist[chain[j]][chain[j+1]]
00089               >dist[chain[i-1]][chain[j]]+dist[chain[i]][chain[j+1]]
00090               +tolerance) {
00091               reverse(i,j); modif=1;
00092            }
00093        }
00094     }
00095        /* Point relink */
00096     for(i=1;i<pointcnt-1;i++) {
00097        double dd;
00098        int th;
00099        th=chain[i];
00100        dd=dist[chain[i-1]][th]+dist[th][chain[i+1]]
00101          -dist[chain[i-1]][chain[i+1]];
00102        for(j=0;j<pointcnt-1;j++) {
00103            if(j==i || j==i-1) continue;
00104            if(dd>dist[chain[j]][th]+dist[chain[j+1]][th]
00105               -dist[chain[j]][chain[j+1]]+tolerance) {
00106               reinsert(i,j+1); modif=1;
00107            }
00108        }
00109     }
00110     for(i=0;i<pointcnt;i++) shpath[shortcnt][i]=chain[i];
00111     shortcnt+=modif; return modif;
00112 }
00113 
00114 void fastcalc(void)
00115 {
00116     int i,j;
00117     
00118     for(i=0;i<pointcnt;i++) chain[i]=i;
00119     shortcnt=1;
00120     for(i=0,j=1;i<10 && j!=0;i++) j=_fcalc();
00121     
00122     for(shortest=0,i=0;i<pointcnt-1;i++) shortest+=dist[chain[i]][chain[i+1]];
00123     if(style==0) shortest+=dist[chain[0]][chain[pointcnt-1]];
00124     for(i=0;i<pointcnt;i++) shpath[0][i]=chain[i];
00125 }
00126 
00127 void _calc(int level, double dis, long used, int utab[])
00128 {
00129     int current[MAX_POINTS];
00130     double dd,d1,d3,d2,dthis;
00131     int i;
00132     long l;
00133 
00134     dd=dis; d1=d2=d3=0;
00135     if(level<pointcnt-2) {
00136        int j,b,prec;
00137        double btab[MAX_POINTS];
00138        if(level>0) prec=utab[level-1]; else prec=0;
00139        if(level>2) {
00140            for(i=0;i<level-2;i++)
00141              btab[i]=dist[utab[i]][prec];
00142            if(level>3) {
00143               d1=ds[level-4]+ds[level-3]+ds[level-2];
00144               d3=dist[utab[level-4]][utab[level-2]]
00145                 +ds[level-2]
00146                 +dist[prec][utab[level-3]];
00147               d2=dist[utab[level-4]][prec]
00148                 +dist[prec][utab[level-3]]
00149                 +ds[level-3];
00150            }
00151        }
00152        if(level==0 || (level==1 && style==0)) b=last+1; else b=0;
00153        for(i=b,l=(long) 1<<b;i<pointcnt;i++,l<<=1) {
00154            if(used&l) continue;
00155            dthis=dist[prec][i];
00156   if(optimize) {
00157            if(style==0) {
00158               if(level>1 && 
00159                  firstlast+dthis>dist[prec][last]+dist[first][i])
00160                 continue;
00161            }
00162            else if(level>0) {
00163                      /* cut differently */
00164               if(style==1 && dthis>firstlast) continue;
00165                      /* switch with first */
00166               if((style&1) && level>1 && dist[first][i]<dthis) continue;
00167                      /* switch with last */
00168               if(style<3 && dist[prec][last]<dthis) continue;
00169            }
00170                      /* already too long */
00171            if(level>0 && shortcnt>0 &&
00172               dthis+dis+distmin*(pointcnt-level-1)>shortest)
00173              continue;
00174                      /* four-config optimize */
00175            for(j=0;j<=level-3;j++) {
00176               if(ds[j]+dthis>btab[j]+dist[utab[j+1]][i])
00177                 goto next;
00178            }
00179                      /* five-config optimize */
00180            if(level>3 && 
00181               (d1+dthis>d2+dist[utab[level-2]][i] ||
00182               d1+dthis>d3+dist[utab[level-3]][i]))
00183              continue;
00184                      /* try prec elsewhere in the chain */
00185            if(level>3) {
00186               double dt;
00187               dt=dist[prec][i]+ds[level-2]-dist[utab[level-2]][i];
00188               for(j=0;j<level-3;j++) {
00189                   if(dt>dist[prec][utab[j]]+dist[prec][utab[j+1]]-ds[j])
00190                     goto next;
00191               }
00192            }
00193   }
00194            memmove(current,utab,level*sizeof(utab[0]));
00195            current[level]=i;
00196            if(level>0) {dd=dis+dthis; ds[level-1]=dthis;}
00197            else {
00198               first=i; firstlast=dist[first][last];
00199            }
00200            _calc(level+1,dd,used|l,current);
00201            next: ;
00202        }
00203        
00204     }
00205     else {
00206        count++;
00207        for(i=0,l=(long) 1;i<pointcnt && used&l; i++,l<<=1);
00208        if(i>=pointcnt) {
00209            printf("ERROR: Internal error.\n"); exit(1);
00210        }
00211        utab[pointcnt-2]=i; utab[pointcnt-1]=last;
00212        dis+=dist[utab[pointcnt-3]][i]+dist[i][last];
00213        if(shortest==0 || dis<shortest-tolerance) {
00214            shortest=dis; shortcnt=0;
00215            memmove(shpath[shortcnt++],utab,MAX_POINTS*sizeof(utab[0]));
00216        }
00217        else if(dis<shortest+tolerance && shortcnt<MAX_SHORT)
00218            memmove(shpath[shortcnt++],utab,MAX_POINTS*sizeof(utab[0]));
00219     }
00220 }
00221 
00222        /* main calculation loops */
00223 void calc(void)
00224 {
00225     int i;
00226     long used;
00227     int utab[MAX_POINTS];
00228 
00229     if(fast) {
00230        fastcalc(); return;
00231     }
00232     count=0; shortest=0; shortcnt=0;
00233     memset(utab,0,sizeof(utab));
00234     switch(style) {
00235        case 0: {     /* loop to the start */
00236            first=0; utab[0]=0;
00237            for(i=1;i<pointcnt;i++) {
00238               last=i; used=(long) 1<<i; used|=(long) 1;
00239               firstlast=dist[first][last];
00240               _calc(1,firstlast,used,utab);
00241            }
00242            break;
00243        }
00244        case 1: {     /* arbitrary open path */
00245            for(i=0;i<pointcnt;i++) {
00246               last=i; used=(long) 1<<i;
00247               _calc(0,0,used,utab);
00248            }
00249            break;
00250        }
00251        case 2: {     /* open path with fixed start */
00252            first=0; utab[0]=0;
00253            for(i=1;i<pointcnt;i++) {
00254               last=i; used=(long) 1|(1<<i);
00255               firstlast=dist[0][last];
00256               _calc(1,0,used,utab);
00257            }
00258            break;
00259        }
00260        case 3: {     /* open path with fixed end */
00261            last=0; used=(long) 1;
00262            _calc(0,0,used,utab);
00263            break;
00264        }
00265        case 4: {     /* open path with fixed start and end */
00266            first=0; utab[0]=0;
00267            last=pointcnt-1; used=(long) 1|(1<<last);
00268            firstlast=dist[first][last];
00269            _calc(1,0,used,utab);
00270            break;
00271        }
00272     }
00273 }
00274 
00275        /* Print the result */
00276 void printresult(void)
00277 {
00278     int i,j;
00279 
00280     printf("%f %d %d\n",shortest,shortcnt,count);
00281     for(j=0;j<shortcnt;j++) {
00282        for(i=0;i<pointcnt;i++) printf("%d ",shpath[j][i]+1);
00283        printf("\n");
00284     }
00285 }
00286 
00287        /* prepare distance table */
00288 void makedist(void)
00289 {
00290     int i,j,k;
00291     double d;
00292 
00293     distmin=-1;
00294     for(i=0;i<pointcnt-1;i++) {
00295        for(j=i+1;j<pointcnt;j++) {
00296            d=0;
00297            for(k=0;k<dimension;k++) {
00298               d+=pow(coord[i][k]-coord[j][k],2);
00299            }
00300            d=sqrt(d); dist[i][j]=dist[j][i]=d;
00301            if(distmin<0 || distmin>d) distmin=d;
00302        }
00303     }
00304     tolerance=distmin*0.000001;
00305 }
00306 
00307 int main(int argc, char *argv[])
00308 {
00309     char *p,*p2,*p3;
00310     int dim;
00311     input=getenv("wims_exec_parm");
00312        /* nothing to do if no problem */
00313     if(input==NULL || *input==0) return 0;
00314 
00315     for(p=input+strlen(input);p>input && isspace(*(p-1));p--);
00316     *p=0; dimension=0;
00317     for(pointcnt=0,p=input;*p!=0 && pointcnt<MAX_POINTS;p=p3) {
00318        p2=strchr(p,'\n');
00319        if(p2==NULL) p3=p+strlen(p);
00320        else {p3=p2+1; *p2=0;}
00321        dim=0;
00322        for(p2=strchr(p,','); p2!=NULL && dim<MAX_DIM;
00323            p=p2+1,dim++,p2=strchr(p,',')) {
00324            *p2=0; coord[pointcnt][dim]=strtod(p,NULL);
00325        }
00326        if(dim<MAX_DIM) coord[pointcnt][dim++]=strtod(p,NULL);
00327        if(dim<2) continue;  /* bad line */
00328        if(dimension==0) dimension=dim;
00329        pointcnt++;
00330     }
00331     if(pointcnt<3) {
00332        printf("ERROR: Not enough points.\n"); exit(1);
00333     }
00334     p=getenv("w_shortpath_style");
00335     if(p!=NULL) style=atoi(p);
00336     if(style<0 || style>4) style=0;
00337     makedist();
00338     optimize=1; calc();printresult();
00339 /*    optimize=0; calc();printresult();
00340 */    return 0;
00341 }
00342