Code covered by the BSD License

Highlights from An Approximation Solution for the Prize-Collecting Steiner Tree Problem

• FindTree(G, vp) G is matrix representing a undirected graph (V,E) have none negtive weighted edges, where the G(i,j) is the cost of edge(i,j)
• PCSTr(G, vp, r) G is matrix representing a undirected graph (V,E) have none negtive weighted edges, where the G(i,j) is the cost of edge(i,j)
• View all files

An Approximation Solution for the Prize-Collecting Steiner Tree Problem

Fangzhou Chen (view profile)

18 Jan 2013 (Updated )

A function giving feasible(not optimal) solutions to the Prize-Collecting Steiner Tree Problem

PCSTr(G, vp, r)
```
%% G is matrix representing a undirected graph (V,E) have none negtive weighted edges, where the G(i,j) is the cost of edge(i,j)
%% Each vertex in V is denoted by the index 1,2,3,... size(V)
%% vp is vertex penalites, a list of int
%% r is the root vertex, the output tree will contain the root
function T = PCSTr(G, vp, r)
F = []; %% A tree
Comp(length(vp)).C = []; %% set of components
Comp(length(vp)).w = 0;
Comp(length(vp)).lambda = 1;
Comp(length(vp)).Id = 0;
RootComp = r; %% the id of root components
d = zeros(length(vp));

nComp = length(vp);  %% number of components

VertexBelong = zeros(1, length(vp)); %% this array is used to denote which component a vertex belongs to. This array is used for
%% eg, VertexBelong(3) = 5 means that the Vertex is in the componet5
%% VertexMark = zeros(1, length(vp));  %% each vertex is labelled with certain component

for i=1: length(vp)
Comp(i).w = 0;
Comp(i).C = [i];
Comp(i).Id = i;
if( i == r )
Comp(i).lambda = 0;
else
Comp(i).lambda = 1;
end
VertexBelong(i) = i; %% initially, each vertex is in its own component
end

next_id = length(vp)+1;  %% next availiable id to assign a component

while (1)
%% while there exists a component in Comp such that its lambda = 1, the while loop continues..
flag = 0;
%% print all comp
%{
for i=1: length(Comp)
fprintf('The comp id is %d with lambda %d.', Comp(i).Id, Comp(i).lambda);
disp(Comp(i).C);
end
%}

for i=1: length(Comp)
if Comp(i).lambda == 1;
flag = 1;
break;
end
end

if (flag==0)
break;
end

%% Find an edge minimize e1
e1 = Inf;
Edge_chosn.i = 0;
Edge_chosn.j = 0;
for i=1:length(G(1,:))

for j=(i+1):length(G(1,:))
%% filter off the entries not representing valid edge cost
if not(G(i,j) >= 0) || VertexBelong(i) == VertexBelong(j)
continue;
end

etmp = (G(i,j) - d(i) - d(j))/( Comp( VertexBelong(i) ).lambda + Comp( VertexBelong(j) ).lambda );
if etmp < e1
e1 = etmp;
Edge_chosn.i = i;
Edge_chosn.j = j;
end

end
end

%% Find an edge minimize e2
e2 = Inf;
Ct_idx = 1;
for i=1:length(Comp)

if (Comp(i).lambda == 1)
etmp = sum( vp( Comp(i).C(1:length(Comp(i).C)) ) );
%{
for k = 1: length(Comp(i).C)
etmp = etmp + vp( Comp(i).C(k) );
end
%}
etmp = etmp - Comp(i).w;

if (etmp < e2)
e2 = etmp;
Ct_idx = i;
end

else
continue;
end
end

e = min(e1,e2);

%% update the w(C) for all component  --- implicitly set yc = yc + lambda(C) for all C
for i=1:length(Comp)
Comp(i).w = Comp(i).w + e*Comp(i).lambda;
end

%% implicitly set ys = 0 for all S belongs to V

for i=1:length(Comp(RootComp).C)
d(  Comp(RootComp).C(i) ) =  d( Comp(RootComp).C(i) ) + e*Comp(RootComp).lambda;
end

if (e == e2)
Comp(Ct_idx).lambda = 0;
%%mark the unlabelled vertices of Ct with label Ct
%{
for i=1:length(Comp(Ct_idx).C)
VertexMark( Comp(Ct_idx).C(i) ) = Comp(Ct_idx).Id;
end
%}
%% disp('not add edge this round');
else
F = [F, Edge_chosn];
Comp(nComp+1).C = [Comp( (VertexBelong(Edge_chosn.i))).C  Comp( (VertexBelong(Edge_chosn.j))).C ];
Comp(nComp+1).w = Comp( (VertexBelong(Edge_chosn.i))).w + Comp( (VertexBelong(Edge_chosn.j))).w;
Comp(nComp+1).Id = next_id;
next_id = next_id + 1;

%%disp(VertexBelong(r));
if( (VertexBelong(r) == (VertexBelong(Edge_chosn.i))) ||  (VertexBelong(r) == (VertexBelong(Edge_chosn.j))) )
Comp(nComp+1).lambda = 0;
else
Comp(nComp+1).lambda = 1;
end

Comp([(VertexBelong(Edge_chosn.i)) (VertexBelong(Edge_chosn.j))]) = [];
nComp = length(Comp);

%% fprintf('add new edge %d %d with lambda = %d', Edge_chosn.i, Edge_chosn.j, Comp(nComp).lambda);
for i=1:length(Comp(nComp).C)
%{
%% update the marks
if VertexMark(Comp(nComp).C(i)) > 0
VertexMark(Comp(nComp).C(i)) = Comp(nComp).Id;
end
%}
VertexBelong(Comp(nComp).C(i)) = nComp;  %% the new added component is at the end of the list
end

%% update the vertex belong
for i=1:length(Comp)
for j=1:length(Comp(i).C)
VertexBelong(Comp(i).C(j)) = i;
if ( Comp(i).C(j) ) == r
RootComp = i;
end
end
end

end

end

% remove the edges!
rm = [];
for i = 1:length(F)
if ( VertexBelong(F(i).i) ~= RootComp && VertexBelong(F(i).j) ~= RootComp )
rm = [rm i];
end
end
F(rm) = [];

%% Generate a structure for the tree!
%% The tree contains two components: 1st is the sets like F; 2nd is the set of vertex objects

vt = [];
T.E = [];
for i = 1:length(F)
vt = union(vt, [F(i).i,F(i).j]);
edge_tmp.id = F(i);
edge_tmp.cost = G(F(i).i, F(i).j);
T.E = [T.E, edge_tmp];
end

%% The vertex objects contain 1. vertex id; 2.vertex neighbours; 3.the revenue of the vertex; 4.whether the vertex is done
T.V = [];

for idx = 1:length(vt)
v_tmp.id = vt(idx);
v_tmp.neighbour = [];
v_tmp.r = vp(v_tmp.id);
v_tmp.done = 0;  %% done is used for pruning the output tree. But currently we do not need it.
%%find the neighbour
for k = 1:length(F)

if(F(k).i == v_tmp.id)
v_tmp.neighbour = union(v_tmp.neighbour, F(k).j);
elseif (F(k).j == v_tmp.id)
v_tmp.neighbour = union(v_tmp.neighbour, F(k).i);
end

end
T.V = [T.V, v_tmp];

end

MapV = zeros(1,length(vp)); %% Map(i) is the index of Vertex i in the tree T.V;  MapV( T.V(i).id ) = i;
for i = 1:length(T.V)
MapV(T.V(i).id) = i;
end

%% After some research, I found that we do not need the strong prune in this example
%% Now we use a recurion method SPT (strong prune tree) to remove as much as edges in T to get the final answers
%{
%% Strong prone tree. But currently we do not need it.
T = SPT(T,G, MapV(T.V(1).id), MapV);

rm = [];
for i = 1: length(T.V)
if T.V(i).done == 3
rm = [rm i];
end
end

T.V(rm) = [];

%}

score = 0;
for i = 1:length(T.V)
T.V(i).r = vp(T.V(i).id);
score = score + vp(T.V(i).id);
end

for i = 1:length(T.E)
score = score - T.E(i).cost;
end

T.score = score;

```