Cycle Detection for Kidney Donor Exchange (netsle02)
/************************************************************************/
/* */
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: netsle02 */
/* TITLE: Cycle Detection for Kidney Donor Exchange (netsle02) */
/* PRODUCT: OR */
/* SYSTEM: ALL */
/* KEYS: OR */
/* PROCS: OPTMODEL, PRINT */
/* DATA: */
/* */
/* SUPPORT: UPDATE: */
/* REF: */
/* MISC: Example 2 from the network solver documentation. */
/* */
/************************************************************************/
/* create random graph on n nodes with arc probability p
and uniform(0,1) weight */
%let n = 100;
%let p = 0.02;
data LinkSetIn;
do from = 0 to &n - 1;
do to = 0 to &n - 1;
if from eq to then continue;
else if ranuni(1) < &p then do;
weight = ranuni(2);
output;
end;
end;
end;
run;
%let max_length = 10;
proc optmodel;
/* declare index sets and parameters, and read data */
set <num,num> ARCS;
num weight {ARCS};
read data LinkSetIn into ARCS=[from to] weight;
set<num,num,num> ID_ORDER_NODE;
/* generate all cycles with 2 <= length <= max_length */
solve with NETWORK /
loglevel = moderate
graph_direction = directed
links = (include=ARCS)
cycle = (mode=all_cycles minlength=2 maxlength=&max_length)
out = (cycles=ID_ORDER_NODE)
;
/* extract <cid,from,to> triples from <cid,order,node> triples */
set <num,num,num> ID_FROM_TO init {};
num last init ., from, to;
for {<cid,order,node> in ID_ORDER_NODE} do;
from = last;
to = node;
last = to;
if order ne 1 then ID_FROM_TO = ID_FROM_TO union {<cid,from,to>};
end;
/* solve set packing problem to find maximum weight node-disjoint union
of short directed cycles */
set CYCLES = setof {<c,i,j> in ID_FROM_TO} c;
set ARCS_c {c in CYCLES} = setof {<(c),i,j> in ID_FROM_TO} <i,j>;
set NODES_c {c in CYCLES} = union {<i,j> in ARCS_c[c]} {i,j};
set NODES = union {c in CYCLES} NODES_c[c];
num cycle_weight {c in CYCLES} = sum {<i,j> in ARCS_c[c]} weight[i,j];
/* UseCycle[c] = 1 if cycle c is used, 0 otherwise */
var UseCycle {CYCLES} binary;
/* declare objective */
max TotalWeight
= sum {c in CYCLES} cycle_weight[c] * UseCycle[c];
/* each node appears in at most one cycle */
con node_packing {i in NODES}:
sum {c in CYCLES: i in NODES_c[c]} UseCycle[c] <= 1;
/* call solver */
solve with milp;
/* output optimal solution */
create data Solution from [c]={c in CYCLES: UseCycle[c].sol > 0.5}
cycle_weight;
quit;
%put &_OROPTMODEL_;
proc print data=Solution noobs label;
run;