Code snippets‎ > ‎

### Binary heap

 12/08/2010: Notice the "arraysize" input variable for dijkstra().  I've found that Mata requires setting the initial minimum array size to a pretty large number (several thousand for a network of only a couple hundred nodes) or else the routine exits with an error.  Setting a large minimum array size eats up a lot of memory and is quite slow.  So I've had to abandon using binary heap. /************************************************************************//* Program: dijkstra().                            *//* Function: Generate distance and path matrices for weighted networks.    *//* Input: Boolean 1 if network directed, 0 if undirected.        *//*   Optional argument to suppress status dots.                *//************************************************************************/transmorphic matrix sgl2::dijkstra(  real scalar z, real scalar arraysize, |string scalar dotsoff){  /* Declare variables.                            */  transmorphic matrix A /* Weighted adjacency list umbrella array.    */  transmorphic matrix T /* Tree array.                    */  transmorphic matrix R /* Return array from bfs().            */  real matrix B /* BFS distance matrix.                    */  real matrix X /* Index matrix.                    */  real matrix N /* Neighbor matrix holding labels and distances.    */  real matrix D /* Distance matrix.                    */  real vector b /* Vector of disconnected nodes.            */  real vector e /* Vector of current estimates.                */  real vector c /* Vector to hold certain/uncertain.            */  real vector v /* Root vector of initial neighbors.            */  real scalar s /* Source vertex.                    */  real scalar m /* Total number of nodes in network.            */  real scalar j /* Column counter.                    */  real scalar exit /* While loop exit indicator.            */  string scalar o /* Status display string.                */  /**********************************************************************/  /* Begin program.                            */  /**********************************************************************/  /* Obtain specified adjacency list and run BFS to obtain information    */  /* on components.                            */  if(z==0){    A=weighted_adjacency_list(0)    R=bfs(0,dotsoff); B=asarray(R,2)  }  else if(z==1){    A=weighted_adjacency_list(1)    R=bfs(1,dotsoff); B=asarray(R,2)  }      /* Obtain total number of nodes in network.                */  m=asarray(A,(3,1))  /* Inform user of program starting.                    */  o="noisily display "+`""""'; stata(o); displayflush()  o="Dijkstra's search algorithm ("+strofreal(m)+" nodes)"  o="noisily display "+`"""'+o+`"""'  stata(o); displayflush()  /* Display dots unless specified not to.                */  if(dotsoff==""){    o="----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5"    o="noisily display "+`"""'+o+`"""'    stata(o); displayflush()   }    /* Define distance matrix.                        */  D=J(m,m,.)  /***** Loop through all nodes.                    */  for(s=1;s<=m;s++){    /* Define vector to hold current estimates.              */                e=J(1,m,.)    /* Define vector to note whether we're certain or not.        */    c=J(1,m,0)    /* Mark distance as being certain.                    */    e[s]=0; c[s]=1       /***** Preliminary check: Continue if source has no neighbors.    */    if(cols(asarray(A,(1,s)))==0){      D[s,.]=e      /* Display status unless specified not to.            */      if(dotsoff=="") display_status(m,s)          continue    }    /***** Mark disconnected components based on info from running BFS.    */    b=select(1..m,B[s,.]:==.)    for(j=1;j<=cols(b);j++) c[b[j]]=1    /***** Generate tree.                        */    /* Define empty tree.                        */    T=asarray_create("real",2,arraysize)    /* Obtain list of neighbors of vertex s sorted by distance.        */    /* No need to add distance - distance of s, or s, is zero.        */    N=(asarray(A,(1,s))\asarray(A,(2,s)))'    _sort(N,2); _transpose(N)           /* Initiate tree.                            */    asarray(T,(1,1),N[.,1])    /* Initiate index matrix: Label, row, row-position.            */    X=(N[1,1],1,1)        /* Loop through each remaining neighbor of s, position them in tree.*/    for(j=2;j<=cols(N);j++){      insert_sift(T,X,N[.,j],s,e)      /* Record neighbors' distance estimates. Since all entries missing*/      /* except for source vertex entry, just straight input.        */      e[N[1,j]]=N[2,j]          }       /* Obtain current root information.                    */    v=asarray(T,(1,1))    /* Mark distance as being certain.                    */    e[v[1]]=v[2]; c[v[1]]=1        /***** Loop while certain/uncertain vector is unfilled.        */    /* Set exit indicator.                        */    exit=0    while(exit!=1){      /***** Exit loop if tree only has a root with no neighbors.    */      if(asarray_elements(T)==1 & cols(asarray(A,(1,v[1])))==0) break      /***** Remove root and sift down if more than one entry in tree.    */      else if(asarray_elements(T)>1) sift_down(T,X)      /***** Add neighbors if v has neighbors.                */      if(cols(asarray(A,(1,v[1])))!=0){    /* Obtain list of neighbors of vertex v.            */    N=(asarray(A,(1,v[1]))\asarray(A,(2,v[1])))         /* Add distance so far from s to v.                */    N[2,.]=N[2,.]:+v[2]        /* Loop through each remaining nbhr of v, position them in tree.*/    for(j=1;j<=cols(N);j++){      insert_sift(T,X,N[.,j],s,e)    }            }       /* Store estimated distance and mark as certain.            */      v=asarray(T,(1,1))      e[v[1]]=v[2]; c[v[1]]=1            /***** Set exit indicator to exit if all distances are certain.    */      if(sum(c)==(cols(c))) exit=1    }      D[s,.]=e    /* Display status unless specified not to.                */    if(dotsoff=="") display_status(m,s)      }  /* Notify program is over.                        */  o="Dijkstra's search algorithm completed"  o="noisily display "+`"""'+o+`"""'  stata(o); displayflush()   /* Return distance matrix.                        */  return(D)}  /************************************************************************//* Program: insert_sift().                        *//* Function: Insert node into next available space in binary tree.    *//* Input: Binary tree, index matrix, 2x1 node vector, and source vertex.*//************************************************************************/void sgl2::insert_sift(transmorphic matrix T, real matrix X,   real vector n, real scalar s, real vector e){  /* Declare variables.                         */  real matrix K /* Array key matrix.                    */  real vector x /* Index input vector.                    */  real vector n2 /* Possible duplication.                 */  real scalar i,j  /* Begin program.                            */  /**********************************************************************/  /* Case: Incoming node is already in tree array.            */  /**********************************************************************/  /* Nodes are in the index matrix only if their estimated distance    */  /* has been recorded in the estimates vector.                */  if(anyof(X[.,1],n[1])){    /* If current estimate is larger than new estimate, update tree,     */    /* index matrix, and estimates vector.                */    if(e[n[1]]>n[2]){      /* First update estimates vector.                    */      e[n[1]]=n[2]      /* Obtain coordinates to current location of node in tree array    */      /* so we can update distance estimates there.            */      n2=select(X,X[.,1]:==n[1])      asarray(T,(n2[2],n2[3]),n)      /* Now that distance estimate has been updated in tree array,     */      /* sift to put node with new estimates in proper place.        */      sift(T,X,n)    }    /* If node is in index matrix (and hence corresponding estimates    */    /* vector is nonmissing) and if new distance is equal to or greater    */    /* than current estimates, don't add new node info and exit function*/    else return  }  /**********************************************************************/  /* Case: Incoming node is not in the tree array.            */  /**********************************************************************/    else if(anyof(X[.,1],n[1])==0){    /* If estimates already defined in the estimates vector, then the    */    /* node has already been in the tree and has been removed with its    */    /* estimate being noted as certain.  Therefore exit function.    */    if(e[n[1]]!=.) return    /* If estimates are not yet defined, record new estimate, and insert*/    /* new node into index matrix and tree array.            */    else if(e[n[1]]==.){      /* Insert estimated distance into estimates vector.        */      e[n[1]]=n[2]      /* Add node into index matrix and tree array.            */      /* Obtain current last index in tree.                */      K=sort(asarray_keys(T),(1,2))      i=K[rows(K),1]; j=K[rows(K),2]      /* Insert node into tree.                        */      if(j+1<=2^(i-1)){    asarray(T,(i,j+1),n)    /* Insert label, row, position.                    */    x=n[1,1]; x=(x,i); x=(x,j+1); X=(X\x)      }         else if(j+1>2^(i-1)){    asarray(T,(i+1,1),n)    /* Insert label, row, position.                    */    x=n[1,1]; x=(x,i+1); x=(x,1); X=(X\x)         }      /***** Sift.                            */      sift(T,X,n)    }    }}/************************************************************************//* Program: sift().                            *//* Function: Sift tree to account for newly added node.            *//* Input: Binary tree and index matrix.                    *//************************************************************************/void sgl2::sift(transmorphic matrix T, real matrix X, real vector n){  /* Declare variables.                         */  pointer(real vector) scalar p1 /* Initial position.            */  pointer(real vector) scalar p2 /* Parent position.             */  real vector x /* Node n's index row.                    */  real vector t /* Parent vector.                    */  real scalar c /* Rounded-up number.                    */  real scalar p /* Position in row.                     */  real scalar r /* Row in tree.                        */  real scalar i1 /* Row index.                        */  real scalar i2 /* Row index.                        */  /* Begin program.                            */  /* Obtain own position.                         */  x=select(X,X[.,1]:==n[1,1])  i1=select(1::rows(X),X[.,1]:==n[1,1])  r=x[1,2]; p=x[1,3]  /* Position pointer at initial start point.                */  p1=&asarray(T,(r,p))    /* Loop while rows are positive.                    */  while(r>1){    /* Find and obtain parent information.                */    c=ceil(p/2)    t=asarray(T,(r-1,c)); p2=&asarray(T,(r-1,c))    /* Compare distance with parent node.                 */    if(t[2]>n[2]){       /* Flip location with parent if distance smaller.            */      *p2=n; *p1=t      /* Update index matrix.                        */      switch_entries(X,n[1],t[1])    }    else if(t[2,1]<=n[2,1]) break    /* Set p1 as current p2.                        */    p1=&asarray(T,(r-1,c))    /* Set new row, position, and index matrix row counter.        */    r=r-1; p=c; i1=i2  }}    /************************************************************************//* Program: sift_down().                        *//* Function: Deletes the smallest value in a tree, brings up the bottom *//*   most entry to the top, then sifts down.                *//* Input: Tree, index matrix, root pointer.                *//************************************************************************/void sgl2::sift_down(transmorphic matrix T, real matrix X){  /* Declare variables.                         */  pointer(real vector) scalar p1 /* Current position.             */  pointer(real vector) scalar left /* Left node pointer.        */  pointer(real vector) scalar right /* Right node pointer.        */   pointer(real vector) scalar root /* Root pointer.            */  real vector r /* Right vector.                    */  real vector l /* Left vector.                        */  real vector n /* Node being moved.                    */  real vector t /* Original root vector.                */  real vector z /* Bottom most vector.                    */  real scalar i /* Row counter.                        */  real scalar p /* Row position.                    */  real scalar x,y /* Coordinate of bottom most node.            */  /* Begin program.                            */  /* Point root pointer to root of tree array.                */  root=&asarray(T,(1,1))  /* Obtain root vector.                        */  t=*root  /* Obtain bottom most vector in tree.                    */  x=max(X[.,2]); y=max(select(X[.,3],X[.,2]:==x))  z=asarray(T,(x,y))  /* Switch entries in index matrix.                    */  switch_entries(X,t[1],z[1])  /* Switch entries in tree array.                    */  *root=asarray(T,(x,y))  /* Remove entry of former root in index matrix.            */  X=select(X,X[.,1]:!=t[1])  /* Remove tree entry of former bottom most node.            */  asarray_remove(T,(x,y))  /***** Begin sifting down.                        */  /* Obtain current node information.                    */  p=1; n=asarray(T,(1,1))  /* Loop until current entry is sifted down into place.        */  for(i=1;i<=(max(X[.,2])-1);i++){    /* Assign left and right pointers.                    */    p1=&asarray(T,(i,p))    left=&asarray(T,(i+1,(p*2)-1))    right=&asarray(T,(i+1,p*2))    /* If both lower right and left nodes are empty, stop sifting down.    */    if(cols(*left)==0 & cols(*right)==0) break    /* If lower left entry is empty but if lower right entry is not    */    /* empty, compare distances.                    */    else if(cols(*left)==0 & cols(*right)!=0){      r=*right      /* Switch tree and index matrix entries if lower entry has    */      /* shorter distance.                         */      if(n[2]>r[2]){        /* Switch tree entries - move bottom entry up then move n down. */        *p1=*right        *right=n        /* Switch index matrix entries.                    */        switch_entries(X,n[1],r[1])        /* Adjust position counter.                    */        p=select(X,X[.,1]:==n[1])[3]              }      else break    }    /* If lower right entry is empty but if lower left entry is not    */    /* empty, compare distances.                    */    else if(cols(*right)==0 & cols(*left)!=0){      l=*left      if(n[2]>l[2]){        /* Switch tree entries - move bottom entry up then move n down. */        *p1=*left        *left=n        /* Switch index matrix entries.                    */        switch_entries(X,n[1],l[1])        /* Adjust position counter.                    */        p=select(X,X[.,1]:==n[1])[3]      }      else break    }        /* If both lower right and left entries are nonempty, compare    */    /* distance with both entries.                    */    else if(cols(*right)!=0 & cols(*left)!=0){      r=*right; l=*left      /* Compare right node.                        */      if(n[2]>r[2] & n[2]<=l[2]){        /* Switch tree entries - move bottom entry up then move n down. */        *p1=*right        *right=n        /* Switch index matrix entries.                    */        switch_entries(X,n[1],r[1])        /* Adjust position counter.                    */        p=select(X,X[.,1]:==n[1])[3]          /* Continue with next iteration.                */        continue      }            /* If distance not larger than right node, compare left node.    */      else if(n[2]>l[2] & n[2]<=r[2]){        /* Switch tree entries - move bottom entry up then move n down. */        *p1=*left        *left=n        /* Switch index matrix entries.                    */        switch_entries(X,n[1],l[1])        /* Adjust position counter.                    */        p=select(X,X[.,1]:==n[1])[3]          /* Continue with next iteration.                */        continue      }      /* If distance is larger than both left and right, switch with    */      /* smaller of the two.                        */      else if(n[2]>l[2] & n[2]>r[2]){        if(l[2]>=r[2]){      /* Switch tree entries - move bottom entry up then move n down*/      *p1=*right      *right=n      /* Switch index matrix entries.                    */      switch_entries(X,n[1],r[1])      /* Adjust position counter.                    */      p=select(X,X[.,1]:==n[1])[3]        /* Continue with next iteration.                */      continue       }    else if(l[2]