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