/* Copyright (C) 1998-2003 XIAO, Gang of Universite de Nice - Sophia Antipolis
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/
/* Finds the shortest paths linking given points. */
/* fast algorithm is not ready. */
#define MAX_POINTS 32 /* cannot exceed number of bits in long */
#define MAX_DIM 10 /* dimension limit */
#define MAX_SHORT 20 /* limit of registered shortest paths */
#include "../config.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <ctype.h>
int pointcnt,dimension,first,last,count,shortcnt;
/* style = 0: loop to the start
* 1: arbitrary open path
* 2: open path with fixed start
* 3: open path with fixed end
* 4: open path with fixed start and end */
int style=0;
int fast=0; /* fast algo if 1, no rigor */
int optimize=0;
char *input;
double dist[MAX_POINTS][MAX_POINTS];
double ds[MAX_POINTS];
double coord[MAX_POINTS][MAX_DIM];
double shortest;
double tolerance; /* to compensate floating rounding error in comparison */
double distmin; /* minimum of distance between two points */
double firstlast; /* distance between first and last */
int shpath[MAX_SHORT][MAX_POINTS];
int chain[MAX_POINTS];
/* reverse a subchain */
void reverse(int i, int j)
{
int k,t;
for(k=0;k<(j-i+1)/2;k++) {
t=chain[i+k]; chain[i+k]=chain[j-k];
chain[j-k]=t;
}
}
/* move point i to position j */
void reinsert(int i, int j)
{
int t;
t=chain[i];
if(j>i) {
j--;
memmove(chain
+i
,chain
+i
+1,(j
-i
)*sizeof(chain
[0]));
chain[j]=t;
}
else {
memmove(chain
+j
+1,chain
+j
,(i
-j
)*sizeof(chain
[0]));
chain[j]=t;
}
}
int _fcalc(void)
{
int modif;
int i,j;
modif=0;
/* four-configuration optimize */
for(i=1;i<pointcnt-2;i++) {
for(j=i+1;j<pointcnt-1;j++) {
if(dist[chain[i-1]][chain[i]]+dist[chain[j]][chain[j+1]]
>dist[chain[i-1]][chain[j]]+dist[chain[i]][chain[j+1]]
+tolerance) {
reverse(i,j); modif=1;
}
}
}
/* Point relink */
for(i=1;i<pointcnt-1;i++) {
double dd;
int th;
th=chain[i];
dd=dist[chain[i-1]][th]+dist[th][chain[i+1]]
-dist[chain[i-1]][chain[i+1]];
for(j=0;j<pointcnt-1;j++) {
if(j==i || j==i-1) continue;
if(dd>dist[chain[j]][th]+dist[chain[j+1]][th]
-dist[chain[j]][chain[j+1]]+tolerance) {
reinsert(i,j+1); modif=1;
}
}
}
for(i=0;i<pointcnt;i++) shpath[shortcnt][i]=chain[i];
shortcnt+=modif; return modif;
}
void fastcalc(void)
{
int i,j;
for(i=0;i<pointcnt;i++) chain[i]=i;
shortcnt=1;
for(i=0,j=1;i<10 && j!=0;i++) j=_fcalc();
for(shortest=0,i=0;i<pointcnt-1;i++) shortest+=dist[chain[i]][chain[i+1]];
if(style==0) shortest+=dist[chain[0]][chain[pointcnt-1]];
for(i=0;i<pointcnt;i++) shpath[0][i]=chain[i];
}
void _calc(int level, double dis, long used, int utab[])
{
int current[MAX_POINTS];
double dd,d1,d3,d2,dthis;
int i;
long l;
dd=dis; d1=d2=d3=0;
if(level<pointcnt-2) {
int j,b,prec;
double btab[MAX_POINTS];
if(level>0) prec=utab[level-1]; else prec=0;
if(level>2) {
for(i=0;i<level-2;i++)
btab[i]=dist[utab[i]][prec];
if(level>3) {
d1=ds[level-4]+ds[level-3]+ds[level-2];
d3=dist[utab[level-4]][utab[level-2]]
+ds[level-2]
+dist[prec][utab[level-3]];
d2=dist[utab[level-4]][prec]
+dist[prec][utab[level-3]]
+ds[level-3];
}
}
if(level==0 || (level==1 && style==0)) b=last+1; else b=0;
for(i=b,l=(long) 1<<b;i<pointcnt;i++,l<<=1) {
if(used&l) continue;
dthis=dist[prec][i];
if(optimize) {
if(style==0) {
if(level>1 &&
firstlast+dthis>dist[prec][last]+dist[first][i])
continue;
}
else if(level>0) {
/* cut differently */
if(style==1 && dthis>firstlast) continue;
/* switch with first */
if((style&1) && level>1 && dist[first][i]<dthis) continue;
/* switch with last */
if(style<3 && dist[prec][last]<dthis) continue;
}
/* already too long */
if(level>0 && shortcnt>0 &&
dthis+dis+distmin*(pointcnt-level-1)>shortest)
continue;
/* four-config optimize */
for(j=0;j<=level-3;j++) {
if(ds[j]+dthis>btab[j]+dist[utab[j+1]][i])
goto next;
}
/* five-config optimize */
if(level>3 &&
(d1+dthis>d2+dist[utab[level-2]][i] ||
d1+dthis>d3+dist[utab[level-3]][i]))
continue;
/* try prec elsewhere in the chain */
if(level>3) {
double dt;
dt=dist[prec][i]+ds[level-2]-dist[utab[level-2]][i];
for(j=0;j<level-3;j++) {
if(dt>dist[prec][utab[j]]+dist[prec][utab[j+1]]-ds[j])
goto next;
}
}
}
memmove(current
,utab
,level
*sizeof(utab
[0]));
current[level]=i;
if(level>0) {dd=dis+dthis; ds[level-1]=dthis;}
else {
first=i; firstlast=dist[first][last];
}
_calc(level+1,dd,used|l,current);
next: ;
}
}
else {
count++;
for(i=0,l=(long) 1;i<pointcnt && used&l; i++,l<<=1);
if(i>=pointcnt) {
}
utab[pointcnt-2]=i; utab[pointcnt-1]=last;
dis+=dist[utab[pointcnt-3]][i]+dist[i][last];
if(shortest==0 || dis<shortest-tolerance) {
shortest=dis; shortcnt=0;
memmove(shpath
[shortcnt
++],utab
,MAX_POINTS
*sizeof(utab
[0]));
}
else if(dis<shortest+tolerance && shortcnt<MAX_SHORT)
memmove(shpath
[shortcnt
++],utab
,MAX_POINTS
*sizeof(utab
[0]));
}
}
/* main calculation loops */
void calc(void)
{
int i;
long used;
int utab[MAX_POINTS];
if(fast) {
fastcalc(); return;
}
count=0; shortest=0; shortcnt=0;
switch(style) {
case 0: { /* loop to the start */
first=0; utab[0]=0;
for(i=1;i<pointcnt;i++) {
last=i; used=(long) 1<<i; used|=(long) 1;
firstlast=dist[first][last];
_calc(1,firstlast,used,utab);
}
break;
}
case 1: { /* arbitrary open path */
for(i=0;i<pointcnt;i++) {
last=i; used=(long) 1<<i;
_calc(0,0,used,utab);
}
break;
}
case 2: { /* open path with fixed start */
first=0; utab[0]=0;
for(i=1;i<pointcnt;i++) {
last=i; used=(long) 1|(1<<i);
firstlast=dist[0][last];
_calc(1,0,used,utab);
}
break;
}
case 3: { /* open path with fixed end */
last=0; used=(long) 1;
_calc(0,0,used,utab);
break;
}
case 4: { /* open path with fixed start and end */
first=0; utab[0]=0;
last=pointcnt-1; used=(long) 1|(1<<last);
firstlast=dist[first][last];
_calc(1,0,used,utab);
break;
}
}
}
/* Print the result */
void printresult(void)
{
int i,j;
printf("%f %d %d\n",shortest
,shortcnt
,count
);
for(j=0;j<shortcnt;j++) {
for(i
=0;i
<pointcnt
;i
++) printf("%d ",shpath
[j
][i
]+1);
}
}
/* prepare distance table */
void makedist(void)
{
int i,j,k;
double d;
distmin=-1;
for(i=0;i<pointcnt-1;i++) {
for(j=i+1;j<pointcnt;j++) {
d=0;
for(k=0;k<dimension;k++) {
d
+=pow(coord
[i
][k
]-coord
[j
][k
],2);
}
d
=sqrt(d
); dist
[i
][j
]=dist
[j
][i
]=d
;
if(distmin<0 || distmin>d) distmin=d;
}
}
tolerance=distmin*0.000001;
}
int main(int argc, char *argv[])
{
char *p,*p2,*p3;
int dim;
input
=getenv("wims_exec_parm");
/* nothing to do if no problem */
if(input==NULL || *input==0) return 0;
*p=0; dimension=0;
for(pointcnt=0,p=input;*p!=0 && pointcnt<MAX_POINTS;p=p3) {
else {p3=p2+1; *p2=0;}
dim=0;
for(p2
=strchr(p
,','); p2
!=NULL
&& dim
<MAX_DIM
;
p
=p2
+1,dim
++,p2
=strchr(p
,',')) {
*p2
=0; coord
[pointcnt
][dim
]=strtod(p
,NULL
);
}
if(dim
<MAX_DIM
) coord
[pointcnt
][dim
++]=strtod(p
,NULL
);
if(dim<2) continue; /* bad line */
if(dimension==0) dimension=dim;
pointcnt++;
}
if(pointcnt<3) {
}
p
=getenv("w_shortpath_style");
if(p
!=NULL
) style
=atoi(p
);
if(style<0 || style>4) style=0;
makedist();
optimize=1; calc();printresult();
/* optimize=0; calc();printresult();
*/ return 0;
}