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
}