geog.c
上传用户:rrhhcc
上传日期:2015-12-11
资源大小:54129k
文件大小:17k
源码类别:

通讯编程

开发平台:

Visual C++

  1. /* $Id: geog.c,v 1.1 1996/10/04 13:36:46 calvert Exp $
  2.  *
  3.  * geog.c -- routines to construct various flavors of semi-geometric graphs
  4.  *           using the Stanford GraphBase tools.
  5.  *
  6.  */
  7. #include <stdio.h>
  8. /* #include <sys/types.h> for NBBY */
  9. /* instead define as foll */
  10. #ifndef NBBY
  11. #define NBBY 8
  12. #endif
  13. #ifndef FBSD
  14. #include <alloca.h>
  15. #endif
  16. #include <assert.h>
  17. #include <string.h> /* for strchr() */
  18. #include "gb_graph.h"
  19. #include "math.h"
  20. #include "geog.h"
  21. #include "limits.h"
  22. #define STACKMAX 0x400 /* 1K bytes max in stack */
  23. /* Convention for trouble codes:
  24.  *  - GEO adds 1-3,
  25.  *  - HIER adds 4-6,
  26.  *  - TS adds 7-9,
  27.  */
  28. #define GEO_TRBL 1
  29. #define HIER_TRBL 4
  30. #define sub u.G /* domain graph vertex expands into */
  31. #ifndef u_char
  32. typedef unsigned char u_char;
  33. #endif 
  34. #ifndef u_long
  35. typedef unsigned long u_long;
  36. #endif
  37. static char geogId[]="$Id: geog.c,v 1.1 1996/10/04 13:36:46 calvert Exp $";
  38. double
  39. distance(Vertex *u, Vertex *v)
  40. {
  41. double dx, dy;
  42. double foo;
  43.     dx = (double) u->xpos - v->xpos;
  44.     dy = (double) u->ypos - v->ypos;
  45.     foo = (dx*dx) + (dy*dy);
  46.     return sqrt(foo);
  47. }
  48. /* Distance between two vertices, rounded to the nearest integer. */
  49. long
  50. idist(Vertex *u, Vertex *v)
  51. {
  52. double dx, dy;
  53. double foo;
  54.     dx = (double) u->xpos - v->xpos;
  55.     dy = (double) u->ypos - v->ypos;
  56.     foo = (dx*dx) + (dy*dy);
  57.     return (long) rint(sqrt(foo));
  58. }
  59. int
  60. printparms(char *buf,geo_parms *pp)
  61. {
  62.         sprintf(buf,"{%ld,%d,%d,%.3f,%.3f,%.3f}",
  63.                 pp->n, pp->scale,pp->method,
  64.                 pp->alpha,pp->beta, pp->gamma);
  65.         /* with Sys V this would be easier */
  66.         return strlen(buf);
  67. }
  68. /* randomize an array of longs, keeping the mean constant */
  69. /* doesn't let anything get below 1, i.e. won't work with negatives */
  70. void
  71. randomize(long* a, long size, long mean, int iters)
  72. {
  73. register i,indx;
  74.    for (i=0; i<size; i++)       /* initialize */
  75.         a[i] = mean;
  76.    if (size < 2) return; /* no point */
  77.    for (i=0; i<iters; i++) {
  78.         indx = gb_next_rand()%size;
  79.         if (a[indx] > 1) {       /* first decrease one */
  80.             a[indx]--;
  81.             a[gb_next_rand()%size]++;   /* then increase one */
  82.         }
  83.    }
  84. }
  85. /*
  86.  * geo(seed, <param structure>)
  87.  *  Creates a gb_graph of n nodes with edges distributed probabilistically.
  88.  *  These graphs are constructed so as to be easy to display graphically.
  89.  *  0. nodes are placed on a grid "scale" units on a side, at most one
  90.  *     node per grid point.
  91.  *     Each node's (integer) position on the grid is kept in the aux fields
  92.  *     x.I and y.I, which are renamed (in geog.h) to xpos and ypos.
  93.  *      If perturb is nonzero, each node's position will also be
  94.  *      offset by a random amount that is uniformly distributed over the
  95.  *      interval [0,resolution/2*scale].
  96.  *  1. Edge is placed between two nodes by a probabilistic method, which
  97.  *     is determined by the "method" parameter.  Edge is placed with
  98.  *     probability p, where p is calculated by one of the methods below,
  99.  *     using:
  100.  *       alpha, beta, gamma: input parameters,
  101.  *       L is scale * sqrt(2.0): the max distance between two points,
  102.  *       d: the Euclidean distance between two nodes.
  103.  *       e: a random number uniformly distributed in [0,L]
  104.  *
  105.  *     Method 1: (Waxman's RG2, with alpha,beta := beta,alpha)
  106.  *          p = alpha * exp(-e/L*beta)
  107.  *     Method 2: (Waxmans's RG1, with alpha,beta := beta, alpha)
  108.  *          p = alpha * exp(-d/L*beta)
  109.  *     Method 3: (Pure random graph)
  110.  *          p = alpha
  111.  *     Method 4: ("EXP" - another distance varying function)
  112.  *     p = alpha * exp(-d/(L-d))
  113.  *     Method 5: (Doar-Leslie, with alpha,beta := beta,alpha, ke=gamma) 
  114.  *     p = (gamma/n) * alpha * exp(-d/(L*beta))
  115.  *     Method 6: (Locality with two regions)
  116.  *          p = alpha     if d <= L*gamma,
  117.  *          p = beta    if d > L*gamma 
  118.  *
  119.  *  2. Constraints (checked upon entry)
  120.  *       0.0 <=  alpha <= 1.0 [alpha is a probability: bad_specs+1]
  121.  *       0.0 <= beta [beta is nonnegative: bad_specs+3]
  122.  *       0.0 <= gamma [gamma is nonnegative: bad_specs+2]
  123.  *       n <  scale*scale  [enough room for nodes: bad_specs+4]
  124.  *
  125.  *  Note carefully the following:
  126.  *   - Output of this routine depends on sequence of values returned by
  127.  *     gb_next_rand().  If the seed parameter is nonzero, it will be
  128.  *     used to seed the random number generator before createing the graph.
  129.  *     if it is zero, whatever random values are next are used.
  130.  *     Thus, successive calls to this function with the same parameters
  131.  *     do not, in general, return the same graph if seed=0.
  132.  *
  133.  *   - The returned graph is NOT guaranteed to be connected.  Checking
  134.  *     connectivity is the responsibility of the caller.
  135.  *
  136.  */
  137. Graph * 
  138. geo(long seed, geo_parms *pp)
  139. {
  140. double p,d,L,radius,r;
  141. u_char *occ;
  142. register int scale;
  143. unsigned nbytes, index, x, y;
  144. u_long units,pertrange;
  145. int mallocd;
  146. register i,j;
  147. Graph *G;
  148. Vertex *up,*vp;
  149. char namestr[128];
  150.     gb_trouble_code = 0;
  151.     if (seed) /* zero seed means "already init'd */
  152. gb_init_rand(seed);
  153.     scale = pp->scale;
  154.     L = sqrt(2.0) * scale;
  155.     gb_trouble_code = 0;
  156.     if ((pp->alpha <= 0.0) || (pp->alpha > 1.0)) gb_trouble_code=bad_specs+1;
  157.     if (pp->gamma < 0.0)  gb_trouble_code=bad_specs+GEO_TRBL;
  158.     if (pp->beta < 0.0) gb_trouble_code=bad_specs+GEO_TRBL;
  159.     if (pp->n > (scale * scale)) gb_trouble_code=bad_specs+GEO_TRBL;
  160.     if (gb_trouble_code)
  161. return NULL;
  162.     radius = pp->gamma * L; /* comes in as a fraction of diagonal */
  163. #define posn(x,y) ((x)*scale)+(y)
  164. #define bitset(v,i) (v[(i)/NBBY]&(1<<((i)%NBBY)))
  165. #define setbit(v,i) v[(i)/NBBY] |= (1<<((i)%NBBY))
  166. /* create a new graph structure */
  167.    
  168.     if ((G=gb_new_graph(pp->n)) == NULL)
  169. return NULL;
  170.     nbytes = ((scale*scale)%NBBY ? (scale*scale/NBBY)+1 : (scale*scale)/NBBY);
  171.     if (nbytes < STACKMAX) { /* small amount - just do it in the stack */
  172.         occ = (u_char *) alloca(nbytes);
  173.         mallocd = 0;
  174.     } else {
  175.         occ = (u_char *) malloc(nbytes);
  176.         mallocd = 1;
  177.     }
  178.     for (i=0; i<nbytes; i++) occ[i] = 0;
  179. /* place each of n points in the plane */
  180.     for (i=0,vp = G->vertices; i<pp->n; i++,vp++) {
  181. sprintf(namestr,"%d",i); /* name is just node number */
  182. vp->name = gb_save_string(namestr);
  183.         do {
  184.             vp->xpos = gb_next_rand()%scale;         /* position: 0..scale-1 */
  185.             vp->ypos = gb_next_rand()%scale;         /* position: 0..scale-1 */
  186.     index = posn(vp->xpos,vp->ypos);
  187.         } while (bitset(occ,index));
  188.         setbit(occ,index);
  189.     }
  190. /* Now go back and add the edges */
  191.     for (i=0,up = G->vertices; i<(pp->n)-1; i++,up++)
  192. for (j=i+1, vp = &G->vertices[i+1]; j<pp->n; j++,vp++) { 
  193.     d = distance(up,vp);
  194.             switch (pp->method) {
  195.     case RG2: /* Waxman's */
  196. d = randfrac()*L;
  197. /* fall through */
  198.             case RG1: /* Waxman's */
  199. p = pp->alpha * exp(-d/(L*pp->beta));
  200. break;
  201.     case RANDOM: /* the classic */
  202. p = pp->alpha;
  203. break;
  204.     case EXP:
  205. p = pp->alpha * exp(-d/(L-d));
  206. break;
  207.             case DL: /* Doar-Leslie */
  208. p = pp->alpha * (pp->gamma/pp->n) * exp(-d/(L*pp->beta));
  209. break;
  210.     case LOC: /* Locality model, two probabilities */
  211. if (d <= radius) p = pp->alpha;
  212. else p = pp->beta;
  213. break;
  214.     default:
  215. die("Bad edge method in geo()!");
  216.     }
  217.     if (randfrac() < p) 
  218. gb_new_edge(up,vp,(int)rint(d));
  219. }
  220. /* Fill in the "id" string to say how the graph was created */
  221.     G->Gscale = scale;
  222.     sprintf(G->id,"geo(%ld,", seed);
  223.     i = strlen(G->id);
  224.     i += printparms(G->id+i,pp);
  225.     G->id[i] = ')';
  226.     G->id[i+1] = (char) 0;
  227.     strcpy(G->util_types,GEO_UTIL);
  228. /* clean up */
  229.     if (mallocd) free(occ);
  230.     return G;
  231. }
  232. #define NOLEAF 0
  233. #define LEAFOK 1
  234. /* Return pointer to node with least degree (least degree > 1, if lflag==0)
  235.  * Such a node is guaranteed to exist, provided the graph is connected and
  236.  * has > 2 nodes.
  237.  */
  238. Vertex *
  239. find_small_deg(Graph *g,int lflag)
  240. {
  241. int i,deg,least=INT_MAX;
  242. Arc *ep;
  243. Vertex *vp,*foo;
  244.    foo = NULL;
  245.    for (i=0,vp=g->vertices; i<g->n; i++,vp++) {
  246. deg = 0;
  247. for (ep=vp->arcs;ep;ep=ep->next) deg++;
  248. if (deg < least) 
  249.     if (lflag == LEAFOK || deg > 1) {
  250. least = deg;
  251. foo = vp;
  252.     }
  253.    }
  254.    return foo;
  255. }
  256. /* find_thresh_deg returns the lowest-numbered node with degree less than
  257.  * threshold, IF such exists.  If all node degrees exceed threshold,
  258.  * returns lowest-numbered node, i.e. g->vertices.
  259.  */
  260. Vertex *
  261. find_thresh_deg(Graph *g,int thresh)
  262. {
  263. int i,deg;
  264. Arc *ep;
  265. Vertex *vp;
  266.    for (i=0,vp=g->vertices; i<g->n; i++,vp++) {
  267. deg = 0;
  268.         for (ep=vp->arcs;ep;ep=ep->next) deg++;
  269. if (deg < thresh)
  270.     return vp;
  271.    }
  272.    return g->vertices;
  273. }
  274. /* Determines what "level" in the graph the edge between u and v is,
  275.  * based on the longest common substring of their "name" strings.  
  276.  */
  277. int
  278. edge_level(Vertex *u, Vertex *v, int nlev)
  279. {
  280. char ss[MAXNAMELEN], tt[MAXNAMELEN];
  281. register char *s=ss, *t=tt, *b, *c;
  282. nlev -= 1;
  283. strcpy(ss,u->name);
  284. strcpy(tt,v->name);
  285. b = strchr(s,'.');
  286. c = strchr(t,'.');
  287. if (b && c) {
  288.   *b++ = '';
  289.   *c++ = '';
  290. }
  291. while (!strcmp(s,t)) { /* still equal => must still be a '.' */
  292.   s = b;
  293.   t = c;
  294.   b = strchr(s,'.');
  295.   c = strchr(t,'.');
  296.   if (b && c) {
  297.     *b++ = '';
  298.     *c++ = '';
  299.   }
  300.   nlev -= 1;
  301. }
  302. return nlev;
  303. /* NOTREACHED */
  304. }
  305. Graph *
  306. geo_hier(long seed,
  307. int nlevels, /* number of levels (=size of following array) */
  308. int edgemeth, /* method of attaching edges */
  309. int aux, /* auxiliary parameter for edge method (threshold) */
  310. geo_parms *pp) /* array of parameter structures, one per level */
  311. {
  312. Graph *newG, *tG, *GG, *srcG, *dstG;
  313. long *numv; /* array of sizes of lower-level graphs */
  314. geo_parms *curparms, workparms[MAXLEVEL];
  315. register i,k,indx;
  316. long dst;
  317. int temp,total,lowsize,otherend,blen,level;
  318. long maxP[MAXLEVEL], maxDiam[MAXLEVEL], wt[MAXLEVEL];
  319. Vertex *np,*vp,*up,*base;
  320. Arc *ap;
  321. char vnamestr[MAXNAMELEN];
  322.      if (seed) /* convention: zero seed means don't use */
  323. gb_init_rand(seed);
  324.      if (nlevels < 1 || nlevels > MAXLEVEL) {
  325. gb_trouble_code = bad_specs+HIER_TRBL;
  326. return NULL;
  327.      }
  328.      /* 1 <= nlevels <= MAXLEVEL */
  329.      /* copy the parameters so we can modify them, and caller doesn't
  330.       * see the changes.
  331.       */
  332.      for (level=0; level<nlevels; level++)
  333. bcopy((char *)&pp[level],&workparms[level],sizeof(geo_parms));
  334.      level = 0;
  335.      gb_trouble_code = 0;
  336.      
  337.      tG = NULL;
  338.      do {
  339. gb_recycle(tG);
  340.         tG = geo(0L,workparms);
  341.      } while (tG != NULL && !isconnected(tG));
  342.      if (tG==NULL)
  343. return tG;
  344.      maxDiam[0] = fdiam(tG);
  345.      maxP[0] = maxDiam[0];
  346.      wt[0] = 1;
  347.      for (i=1; i<nlevels; i++)
  348.        maxDiam[i] = -1;
  349.      curparms = workparms;
  350.      while (++level < nlevels) {
  351.        long tdiam;
  352. curparms++; /* parameters for graphs @ next level */
  353.         /* spread out the numbers of nodes per graph at this level */
  354. numv = (long *) calloc(tG->n,sizeof(long));
  355. lowsize = curparms->n;
  356.         randomize(numv,tG->n,curparms->n,3*tG->n);
  357. /* create a subordinate graph for each vertex in the "top" graph,
  358.          * and add it into the new graph as a whole.
  359.  * We construct the subgraphs all at once to ensure that each
  360.  * has a unique address.
  361.  */
  362. for (i=0,vp=tG->vertices; i<tG->n; i++,vp++) {
  363.     curparms->n = numv[i];
  364.     do {
  365. newG = geo(0L,curparms);
  366. if (newG==NULL) return NULL;
  367.     } while (!isconnected(newG));
  368.     vp->sub = newG;
  369.     tdiam = fdiam(newG);
  370.     if (tdiam>maxDiam[level])
  371.       maxDiam[level] = tdiam;
  372. }
  373. /* do some calculations before "flattening" the top Graph */
  374. total = 0;
  375. for (i=0; i<tG->n; i++) { /* translate node numbers */
  376.     temp = numv[i];
  377.     numv[i]= total;
  378.     total += temp;
  379. }
  380. if (total != tG->n*lowsize) {
  381.     fprintf(stderr,"bad size of new graph!n");
  382.     fprintf(stderr,"total %d tG->n %d lowsize %dn",total,tG->n,lowsize);
  383.     gb_trouble_code = impossible+HIER_TRBL;
  384.     return NULL;
  385. }
  386. /* now create what will become the "new" top-level graph */
  387. newG = gb_new_graph(total);
  388. if (newG==NULL) {
  389.     gb_trouble_code += HIER_TRBL;
  390.     return NULL;
  391. }
  392. /* resolution of the new graph */
  393. newG->Gscale = tG->Gscale * curparms->scale;
  394.        /* compute edge weights for this level */
  395.        wt[level] = maxP[level-1] + 1;
  396.        maxP[level] = (maxDiam[level]*wt[level])
  397.  + (maxDiam[level-1]*maxP[level-1]);
  398.        for (i=0,vp=tG->vertices; i<tG->n; i++,vp++) {
  399.  strcpy(vnamestr,vp->name); /* base name for all "offspring" */
  400.  blen = strlen(vnamestr);
  401.  vnamestr[blen] = '.';
  402.     
  403.  GG = tG->vertices[i].sub;
  404.  base = newG->vertices + numv[i]; /* start of this node's */
  405.  for (k=0,np=base,up=GG->vertices; k<GG->n; k++,np++,up++) {
  406.    /* add the node's edges */
  407.    for (ap=up->arcs; ap; ap=ap->next)  {
  408.      otherend = ap->tip - GG->vertices;
  409.      if (k < otherend)
  410.        gb_new_edge(np,base+otherend,ap->len);
  411.      
  412.    }
  413.    /* now set the new node's position */
  414.    np->xpos = tG->vertices[i].xpos * curparms->scale + up->xpos;
  415.    np->ypos = tG->vertices[i].ypos * curparms->scale + up->ypos;
  416.    /* give the "new" node a name by catenating top & bot names */
  417.    strcpy(vnamestr+blen+1,up->name);
  418.    np->name = gb_save_string(vnamestr);
  419.  } /* loop over GG's vertices */
  420.        }  /* loop over top-level vertices */
  421.        /*
  422. * Now we have to transfer the top-level edges to new graph.
  423. * This is done by one of three methods:
  424. *    0: choose a random node in each subgraph
  425. *    1: attach to the smallest-degree non-leaf node in each
  426. *    2: attach to smallest-degree node
  427. *    3: attach to first node with degree less than aux
  428. */
  429.        for (i=0; i<tG->n; i++) {
  430.  Vertex *srcp, *dstp;
  431.  Graph *srcG, *dstG;
  432.  srcG = tG->vertices[i].sub;
  433.  if (srcG == NULL) { /* paranoia */
  434.    gb_trouble_code = impossible+HIER_TRBL+1;
  435.    return NULL;
  436.  }
  437.  for (ap=tG->vertices[i].arcs; ap; ap=ap->next) {
  438.    dst = ap->tip - tG->vertices;
  439.    if (i > dst) /* consider each edge only ONCE */
  440.      continue;
  441.    dstG = ap->tip->sub;
  442.    if (dstG == NULL) { /* paranoia */
  443.      gb_trouble_code = impossible+HIER_TRBL+1;
  444.      return NULL;
  445.    }
  446.    /* choose endpoints of the top-level edge */
  447.    switch (edgemeth) {
  448.    case 0: /* choose random node in each */
  449.      srcp = srcG->vertices + gb_next_rand()%srcG->n;
  450.      dstp = dstG->vertices + gb_next_rand()%dstG->n;
  451.      break;
  452.    case 1: /* find nonleaf node of least degree in each */
  453.      /* This causes problems with graph size < 3 */
  454.      if (srcG->n > 2)
  455.        srcp = find_small_deg(srcG,NOLEAF);
  456.      else
  457.        srcp = find_small_deg(srcG,LEAFOK);
  458.      if (dstG->n > 2)
  459.        dstp = find_small_deg(dstG,NOLEAF);
  460.      else
  461.        dstp = find_small_deg(dstG,LEAFOK);
  462.      break;
  463.    case 2: /* find node of smallest degree */
  464.      srcp = find_small_deg(srcG,LEAFOK);
  465.      dstp = find_small_deg(dstG,LEAFOK);
  466.      break;
  467.    case 3: /* first node w/degree < aux */
  468.      srcp = find_thresh_deg(srcG,aux);
  469.      dstp = find_thresh_deg(dstG,aux);
  470.    default:
  471.      gb_trouble_code = bad_specs+HIER_TRBL;
  472.      return NULL;
  473.    } /* switch on edgemeth */
  474.    /* pointer arithmetic: isn't it fun?
  475.       printf("Copying edge from %d to %dn",
  476.       numv[i]+(srcp - srcG->vertices),
  477.       numv[dst] + (dstp - dstG->vertices));
  478.       */
  479.    if (srcp==NULL || dstp==NULL) {
  480.      gb_trouble_code = impossible + HIER_TRBL+2;
  481.      return NULL;
  482.    }
  483.    srcp = newG->vertices + numv[i] + (srcp - srcG->vertices);
  484.    dstp = newG->vertices + numv[dst] + (dstp - dstG->vertices);
  485.    gb_new_edge(srcp,dstp,idist(srcp,dstp));
  486.  } /* for each arc */
  487.        } /* for each vertex of top graph */
  488.         /* now make the "new" graph the "top" graph and recycle others */
  489.        for (i=0,vp=tG->vertices; i<tG->n; i++,vp++)
  490.  gb_recycle(vp->sub);
  491.        gb_recycle(tG);
  492.        tG = newG;
  493.        free(numv);
  494.      } /* while more levels */
  495. /* Finally, go back and add the policy weights,
  496.  * based upon the computed max diameters
  497.  * and Max Path lengths.
  498.  */
  499.    for (i=0; i<tG->n; i++)
  500.      for (ap=tG->vertices[i].arcs; ap; ap=ap->next) {
  501.        dst = ap->tip - tG->vertices;
  502.        if (i > dst) /* consider each edge only ONCE */
  503.  continue;
  504.        assert(i != dst); /* no self loops */
  505.        /* i < dst: it is safe to refer to ap's mate by ap+1.  */
  506.        level = edge_level(&tG->vertices[i],&tG->vertices[dst],nlevels);
  507.        ap->policywt = (ap+1)->policywt = wt[level];
  508.      }
  509. /* construct the utility and id strings for the new graph.
  510.  * Space constraints will restrict us to keeping about 4 levels'
  511.  * worth of info.
  512.  */
  513.   {
  514.     char buf[ID_FIELD_SIZE+1];
  515.     register char *cp;
  516.     int len, nextlen, left;
  517.     strcpy(tG->util_types,GEO_UTIL); /* same for all geo graphs, */
  518. /* defined in geo.h */
  519.     cp = tG->id;
  520.     sprintf(cp,"geo_hier(%d,%d,%d,%d,[",seed,nlevels,edgemeth,aux);
  521.     len = strlen(cp);
  522.     left = ID_FIELD_SIZE - len;
  523.     cp += len;
  524.     
  525.     for (i=0; (i < nlevels) && (left > 0); i++) {
  526.       nextlen = printparms(buf,&pp[i]);
  527.       strncpy(cp,buf,left);
  528.       left -= nextlen;
  529.       cp += nextlen;
  530.     }
  531.     if (left > 0) {
  532.       sprintf(buf,"])");
  533.       nextlen = strlen(buf);
  534.       strncpy(cp,buf,left);
  535.     }
  536.   }
  537.   return tG;
  538. } /* geo_hier() */