Cycle Detection for Kidney Donor Exchange
/***************************************************************/
/* */
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: onete02 */
/* TITLE: Cycle Detection for Kidney Donor Exchange */
/* PRODUCT: OR */
/* SYSTEM: ALL */
/* KEYS: OR */
/* PROCS: OPTNET, OPTMODEL, PRINT */
/* DATA: */
/* */
/* SUPPORT: UPDATE: */
/* REF: */
/* MISC: Example 2 from the OPTNET 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;
/* generate all cycles with 2 <= length <= max_length */
%let max_length = 10;
proc optnet
loglevel = moderate
graph_direction = directed
data_links = LinkSetIn;
cycle
minLength = 2
maxLength = &max_length
out = Cycles
mode = all_cycles;
run;
%put &_OROPTNET_;
%put &_OROPTNET_CYCLE_;
/* convert Cycles into one observation per arc */
data Cycles0(keep=c i j);
set Cycles;
retain last;
c = cycle;
i = last;
j = node;
last = j;
if order ne 1 then output;
run;
/* solve set packing problem to find maximum weight node-disjoint union
of short directed cycles */
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> TRIPLES;
read data Cycles0 into TRIPLES=[c i j];
set CYCLES = setof {<c,i,j> in TRIPLES} c;
set ARCS_c {c in CYCLES} = setof {<(c),i,j> in TRIPLES} <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;