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]<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         
        }
      }
      /* Distance is < | = to either lower right or left and sift    */
      /* down should be stopped.                    */
      else break
    }   
  }
}
/************************************************************************/
/* Program: switch_entries().                        */
/* Function: Switch row entries of an index matrix.            */
/* Input: Index matrix, label a and label b.                */
/************************************************************************/
void sgl2::switch_entries(real matrix X, real scalar a, real scalar b){
  /* Declare variables.                         */
  real scalar i1 /* Row index of node a.                */
  real scalar i2 /* Row index of node b.                */
  /* Begin program.                             */
  /* Obtain row index of rows corresponding to nodes a and b.        */
  i1=select((1::rows(X)),X[.,1]:==a)
  i2=select((1::rows(X)),X[.,1]:==b)
  /* Switch rows.                            */
  X[i1,1]=b; X[i2,1]=a
}
Comments