wolffd@0: function A = minimum_spanning_tree(C1, C2)
wolffd@0: %
wolffd@0: % Find the minimum spanning tree using Prim's algorithm.
wolffd@0: % C1(i,j) is the primary cost of connecting i to j.
wolffd@0: % C2(i,j) is the (optional) secondary cost of connecting i to j, used to break ties.
wolffd@0: % We assume that absent edges have 0 cost.
wolffd@0: % To find the maximum spanning tree, used -1*C.
wolffd@0: % See Aho, Hopcroft & Ullman 1983, "Data structures and algorithms", p 237.
wolffd@0: 
wolffd@0: % Prim's is O(V^2). Kruskal's algorithm is O(E log E) and hence is more efficient
wolffd@0: % for sparse graphs, but is implemented in terms of a priority queue.
wolffd@0: 
wolffd@0: % We partition the nodes into those in U and those not in U.
wolffd@0: % closest(i) is the vertex in U that is closest to i in V-U.
wolffd@0: % lowcost(i) is the cost of the edge (i, closest(i)), or infinity is i has been used.
wolffd@0: % In Aho, they say C(i,j) should be "some appropriate large value" if the edge is missing.
wolffd@0: % We set it to infinity.
wolffd@0: % However, since lowcost is initialized from C, we must distinguish absent edges from used nodes.
wolffd@0: 
wolffd@0: n = length(C1);
wolffd@0: if nargin==1, C2 = zeros(n); end
wolffd@0: A = zeros(n);
wolffd@0: 
wolffd@0: closest = ones(1,n);
wolffd@0: used = zeros(1,n); % contains the members of U
wolffd@0: used(1) = 1; % start with node 1
wolffd@0: C1(find(C1==0))=inf;
wolffd@0: C2(find(C2==0))=inf;
wolffd@0: lowcost1 = C1(1,:);
wolffd@0: lowcost2 = C2(1,:);
wolffd@0: 
wolffd@0: for i=2:n
wolffd@0:   ks = find(lowcost1==min(lowcost1));
wolffd@0:   k = ks(argmin(lowcost2(ks)));
wolffd@0:   A(k, closest(k)) = 1;
wolffd@0:   A(closest(k), k) = 1;
wolffd@0:   lowcost1(k) = inf;
wolffd@0:   lowcost2(k) = inf;
wolffd@0:   used(k) = 1;
wolffd@0:   NU = find(used==0);
wolffd@0:   for ji=1:length(NU)
wolffd@0:     for j=NU(ji)
wolffd@0:       if C1(k,j) < lowcost1(j)
wolffd@0: 	lowcost1(j) = C1(k,j);
wolffd@0: 	lowcost2(j) = C2(k,j);
wolffd@0: 	closest(j) = k;
wolffd@0:       end
wolffd@0:     end
wolffd@0:   end
wolffd@0: end
wolffd@0: