gb_econ.w

模拟器提供了一个简单易用的平台

strcpy((p+1)->title,"Users");@+node_index[MAX_N]=p+1;

@ The remaining part of \.{econ.dat} is an $81\times80$ matrix in which
the $k$th row contains the outputs of sector~$k$ to all sectors except
\.{Users}. Each row consists of a blank line followed by 8 data lines;
each data line contains 10 numbers separated by commas.
Zeroes are represented by |""| instead of by |"0"|.
For example, the data line $$\hbox{\tt
8490,2182,42,467,,,,,,}$$ follows the initial blank line; it means
that sector~1 output 8490 million dollars to itself, \$2182M to
sector~2, \dots, \$0M to sector~10.

@<Read and store the output...@>=
{@+register long s=0; /* row sum */
  register long x; /* entry read from \.{econ.dat} */
  if (gb_char()!='\n') panic(syntax_error+5);
   /* blank line missing between rows */
  gb_newline();
  p=node_index[k];
  for (j=1;j<MAX_N;j++) {
    p->table[j]=x=gb_number(10);@+s+=x;
    node_index[j]->total+=x;
    if ((j%10)==0) {
      if (gb_char()!='\n') panic(syntax_error+6);
       /* out of synch in input file */
      gb_newline();
    }@+else if (gb_char()!=',') panic(syntax_error+7);
     /* missing comma after entry */
  }
  p->table[MAX_N]=s; /* sum of |table[1]| through |table[80]| */
}

@* Growing a subtree.
Once all the data appears in |node_block|, we want to extract from it
and combine~it as specified by parameters |n|, |omit|, and |seed|.
This amalgamation process effectively prunes the tree; it can also be
regarded as a procedure that grows a subtree of the full economic tree.

@<Determine the |n| sectors to use in the graph@>=
{@+long l=n+omit-2; /* the number of leaves in the desired subtree */
  if (l==NORM_N) @<Choose all sectors@>@;
  else if (seed) @<Grow a random subtree with |l| leaves@>@;
  else @<Grow a subtree with |l| leaves by subdividing largest sectors first@>;
}

@ The chosen leaves of our subtree are identified by having their
|tag| field set to~1.

@<Choose all sectors@>=
for (k=NORM_N;k;k--) node_index[k]->tag=1;

@ To grow the |l|-leaf subtree when |seed=0|, we first pass over the
tree bottom-up to compute the total input (and output) of each macro-sector;
then we proceed from the top down to subdivide sectors in decreasing
order of their total input. This process provides a good introduction to the
bottom-up and top-down tree methods we will be using in several other
parts of the program.

The |special| node is used here for two purposes: It is the head of a
linked list of unexplored nodes, sorted by decreasing order of
their |total| fields; and it appears at the end of that list, because
|special->total=0|.

@<Grow a subtree with |l| leaves by subdividing largest sectors first@>=
{@+register node *special=node_index[MAX_N];
     /* the \.{Users} node at the end of |node_block| */
  for (p=node_index[ADJ_SEC]-1;p>=node_block;p--) /* bottom up */
    if (p->rchild)
      p->total=(p+1)->total+p->rchild->total;
  special->link=node_block;@+node_block->link=special; /* start at the root */
  k=1; /* |k| is the number of nodes we have tagged or put onto the list */
  while (k<l) @<If the first node on the list is a leaf, delete it and tag it;
                otherwise replace it by its two children@>;
  for (p=special->link;p!=special;p=p->link)
    p->tag=1; /* tag everything on the list */
}

@ @<If the first node on the list is a leaf,...@>=
{
  p=special->link; /* remove |p|, the node with greatest |total| */
  special->link=p->link;
  if (p->rchild==0) p->tag=1; /* |p| is a leaf */
  else {
    pl=p+1;@+pr=p->rchild;
    for (q=special;q->link->total>pl->total;q=q->link) ;
    pl->link=q->link;@+q->link=pl; /* insert left child in its proper place */
    for (q=special;q->link->total>pr->total;q=q->link) ;
    pr->link=q->link;@+q->link=pr; /* insert right child in its proper place */
    k++;
  }
}

@ We can obtain a uniformly distributed |l|-leaf subtree of a given tree
by choosing the root when |l=1| or by using the following idea when |l>1|:
Suppose the given tree~$T$ has subtrees $T_0$ and $T_1$. Then it has
$T(l)$ subtrees with |l|~leaves, where $T(l)=\sum_k T_0(k)T_1(l-k)$.
We choose a random number $r$ between 0 and $T(l)-1$, and we find the
smallest $m$ such that $\sum_{k\le m}T_0(k)T_1(l-k)>r$. Then we
proceed recursively to
compute a random $m$-leaf subtree of~$T_0$ and a random $(l-m)$-leaf
subtree of~$T_1$.

A difficulty arises when $T(l)$ is $2^{31}$ or more. But then we can replace
$T_0(k)$ and $T_1(l-k)$ in the formulas above by $\lceil T_0(k)/d_0\rceil$
and $\lceil T_1(k)/d_1\rceil$, respectively, where $d_0$ and $d_1$ are
arbitrary constants; this yields smaller values
$T(l)$ that define approximately the same distribution of~$k$.

The program here computes the $T(l)$ values bottom-up, then grows a
random tree top-down. If node~|p| is not a leaf, its |table[0]| field
will be set to the number of leaves below it; and its |table[l]| field
will be set to $T(l)$, for |1<=l<=table[0]|.

The data in \.{econ.dat} is sufficiently simple that most of the $T(l)$
values are less than $2^{31}$. We need to scale them
down to avoid overflow only at the root node of the tree; this
case is handled separately.

We set the |tag| field of a node equal to the number of leaves to be
grown in the subtree rooted at that node. This convention is consistent
with our previous stipulation that |tag=1| should characterize the
nodes that are chosen to be vertices.

@<Grow a random subtree with |l| leaves@>=
{
  node_block->tag=l;
  for (p=node_index[ADJ_SEC]-1;p>node_block;p--) /* bottom up, except root */
    if (p->rchild) @<Compute the $T(l)$ values for subtree |p|@>;
  for (p=node_block;p<node_index[ADJ_SEC];p++) /* top down, from root */
    if (p->tag>1) {
      l=p->tag;
      pl=p+1;@+pr=p->rchild;
      if (pl->rchild==NULL) {
        pl->tag=1;@+pr->tag=l-1;
      }@+else if (pr->rchild==NULL) {
        pl->tag=l-1;@+pr->tag=1;
      }@+else @<Stochastically determine the number of leaves to grow in
                    each of |p|'s children@>;
    }
}

@ Here we are essentially multiplying two generating functions.
Suppose $f(z)=\sum_l T(l)z^l$; then we are computing $f_p(z)=
z+f_{pl}(z)f_{pr}(z)$.

@<Compute the $T(l)$ values for subtree |p|@>=
{
  pl=p+1;@+pr=p->rchild;
  p->table[1]=p->table[2]=1; /* $T(1)$ and $T(2)$ are always 1 */
  if (pl->rchild==0) { /* left child is a leaf */
    if (pr->rchild==0) p->table[0]=2; /* and so is the right child */
    else { /* no, it isn't */
      for (k=2;k<=pr->table[0];k++) p->table[1+k]=pr->table[k];
      p->table[0]=pr->table[0]+1;
    }
  }@+else if (pr->rchild==0) { /* right child is a leaf */
    for (k=2;k<=pl->table[0];k++) p->table[1+k]=pl->table[k];
    p->table[0]=pl->table[0]+1;
  }@+else { /* neither child is a leaf */
    @<Set |p->table[2]|, |p->table[3]|, \dots\ to convolution of
      |pl| and |pr| table entries@>;
    p->table[0]=pl->table[0]+pr->table[0];
  }
}

@ @<Set |p->table[2]|, |p->table[3]|, \dots\ to convolution...@>=
p->table[2]=0;
for (j=pl->table[0];j;j--) {@+register long t=pl->table[j];
  for (k=pr->table[0];k;k--)
    p->table[j+k]+=t*pr->table[k];
}

@ @<Stochastically determine the number of leaves to grow...@>=
{@+register long ss,rr;
  j=0; /* we will set |j=1| if scaling is necessary at the root */
  if (p==node_block) {
    ss=0;
    if (l>29 && l<67) {
      j=1; /* more than $2^{31}$ possibilities exist */
      for (k=(l>pr->table[0]? l-pr->table[0]: 1);k<=pl->table[0] && k<l;k++)
        ss+=((pl->table[k]+0x3ff)>>10)*pr->table[l-k];
              /* scale with $d_0=1024$, $d_1=1$ */
    }@+else
      for (k=(l>pr->table[0]? l-pr->table[0]: 1);k<=pl->table[0] && k<l;k++)
        ss+=pl->table[k]*pr->table[l-k];
  }@+else ss=p->table[l];
  rr=gb_unif_rand(ss);
  if (j)
    for (ss=0,k=(l>pr->table[0]? l-pr->table[0]: 1);ss<=rr;k++)
      ss+=((pl->table[k]+0x3ff)>>10)*pr->table[l-k];
  else for (ss=0,k=(l>pr->table[0]? l-pr->table[0]: 1);ss<=rr;k++)
      ss+=pl->table[k]*pr->table[l-k];
  pl->tag=k-1;@+pr->tag=l-k+1;
}

@* Arcs.
In the general case, we have to combine some of the basic micro-sectors
into macro-sectors by adding together the appropriate input/output
coefficients. This is a bottom-up pruning process.

Suppose |p| is being formed as the union of |pl| and~|pr|.
Then the arcs leading out of |p| are obtained by summing the numbers
on arcs leading out of |pl| and~|pr|; the arcs leading into |p| are
obtained by summing the numbers on arcs leading into |pl| and~|pr|;
the arcs from |p| to itself are obtained by summing the four numbers
on arcs leading from |pl| or~|pr| to |pl| or~|pr|.

We maintain the |node_index| table so that its non-|NULL| entries
contain all the currently active nodes. When |pl| and~|pr| are
being pruned in favor of~|p|, node |p|~inherits |pl|'s place in
|node_index|; |pr|'s former place becomes~|NULL|.

@<Put the appropriate arcs into the graph@>=
@<Prune the sectors that are used in macro-sectors, and form
  the lists of SIC sector codes@>;
@<Make the special nodes invisible if they are omitted, visible otherwise@>;
@<Compute individual thresholds for each chosen sector@>;
{@+register Vertex *v=new_graph->vertices+n;
  for (k=MAX_N;k;k--)
    if ((p=node_index[k])!=NULL) {
      vert_index[k]=--v;
      v->name=gb_save_string(p->title);
      v->SIC_codes=p->SIC_list;
      v->sector_total=p->total;
    }@+else vert_index[k]=NULL;
  if (v!=new_graph->vertices)
    panic(impossible); /* bug in algorithm; this can't happen */
  for (j=MAX_N;j;j--)
    if ((p=node_index[j])!=NULL) {@+register Vertex *u=vert_index[j];
      for (k=MAX_N;k;k--)
        if ((v=vert_index[k])!=NULL)
          if (p->table[k]!=0 && p->table[k]>node_index[k]->thresh) {
            gb_new_arc(u,v,1L);
            u->arcs->flow=p->table[k];
          }
      }
}

@ @<Private v...@>=
static Vertex *vert_index[MAX_N+1]; /* the vertex assigned to an SIC code */

@ The theory underlying this step is the following, for integers
$a,b,c,d$ with $b,d>0$:
$$ {a\over b}>{c\over d} \qquad\iff\qquad
  a>\biggl\lfloor{b\over d}\biggr\rfloor\,c +
       \biggl\lfloor{(b\bmod d)c\over d}\biggr\rfloor\,.$$
In our case, $b=\hbox{|p->total|}$ and $c=threshold\le d=65536=2^{16}$, hence
the multiplications cannot overflow. (But they can come awfully darn close.)

@<Compute individual thresholds for each chosen sector@>=
for (k=MAX_N;k;k--)
  if ((p=node_index[k])!=NULL) {
    if (threshold==0) p->thresh=-99999999;
    else p->thresh=((p->total>>16)*threshold)+
      (((p->total&0xffff)*threshold)>>16);
  }

@ @<Prune the sectors that are used in macro-sectors, and form
  the lists of SIC sector codes@>=
for (p=node_index[ADJ_SEC];p>=node_block;p--) { /* bottom up */
  if (p->SIC) { /* original leaf */
    p->SIC_list=gb_virgin_arc();
    p->SIC_list->len=p->SIC;
  }@+else {
    pl=p+1;@+pr=p->rchild;
    if (p->tag==0) p->tag=pl->tag+pr->tag;
    if (p->tag<=1) @<Replace |pl| and |pr| by their union, |p|@>;
  }
}
    
@ @<Replace |pl| and |pr| by their union, |p|@>=
{@+register Arc *a=pl->SIC_list;
  register long jj=pl->SIC, kk=pr->SIC;
  p->SIC_list=a;
  while (a->next) a=a->next;
  a->next=pr->SIC_list;
  for (k=MAX_N;k;k--)
    if ((q=node_index[k])!=NULL) {
      if (q!=pl && q!=pr) q->table[jj]+=q->table[kk];
      p->table[k]=pl->table[k]+pr->table[k];
    }
  p->total=pl->total+pr->total;
  p->SIC=jj;
  p->table[jj]+=p->table[kk];
  node_index[jj]=p;
  node_index[kk]=NULL;
}

@ If the \.{Users} vertex is not omitted, we need to compute each
sector's total final demand, which is calculated so that the row sums
and column sums of the input/output coefficients come out equal. We've
already computed the column sum, |p->total|; we've also computed
|p->table[1]+@[@t\hbox{$\cdots$}@>@]+p->table[ADJ_SEC]|, and put it into
|p->table[MAX_N]|. So now we want to replace |p->table[MAX_N]| by
|p->total-p->table[MAX_N]|. As remarked earlier, this quantity might
be negative.

In the special node |p| for the \.{Users} vertex, the preliminary
processing has made |p->total=0|; moreover, |p->table[MAX_N]| is the
sum of value added, or GNP.  We want to switch those fields.

We don't have to set the |tag| fields to 1 in the special nodes, because
the remaining parts of the arc-generation algorithm don't look at those fields.

@<Make the special nodes invisible if they are omitted, visible otherwise@>=
if (omit==2) node_index[ADJ_SEC]=node_index[MAX_N]=NULL;
else if (omit==1) node_index[MAX_N]=NULL;
else {
  for (k=ADJ_SEC;k;k--)
    if ((p=node_index[k])!=NULL) p->table[MAX_N]=p->total-p->table[MAX_N];
  p=node_index[MAX_N]; /* the special node */
  p->total=p->table[MAX_N];
  p->table[MAX_N]=0;
}
      
@* Index. As usual, we close with an index that
shows where the identifiers of \\{gb\_econ} are defined and used.