How to use this notebook

This notebook is meant to solve maximum stable set instances by using polyhedral combinatorics techniques.
The simplex method is applied to the dual problem of the stable set inequality description.
First only arc inequalities are considered  (i.e.  x_i+x_j <= 1) and nonnegativity
In Phase I an initial feasible basis solution is found for its dual. 
Phase II reoptimizes the dual after any added primal inequality (dual column).

The solution is displayed at each iteration.  Nodes  with value 1  are green. Nodes  with value 0  are blue. Nodes with fractional values are red with node size proportional to the value.

Clique inequalities must be found by inspection.

Odd hole inequalities may be found either by manual inspection or automatically.

The process has to be repeated until an integer solution is found or no evident cuts are available. In the latter case branching should be required. At the moment this feature is not implemented.

To use the noteook one has first to run the section Routines. Instances are input  in  the section Data.
There are several ways to input instances: these are described within the section. 
The computation is carried out in the section Phase I. and then in the section  Phase II. This has to be repeated as many times as necessary. Instructions are found within the section.

Routines

Off[General::spell1];

Off[General::spell];

findoddwheel[center_,z_]:=Module[{dist,result,auxarcs,currnodes,arcs,kk1,kk2},
dist={};violations={};
auxarcs={};
Do[{kk1,kk2}=subarcs[[center,i]];
     auxarcs=Join[auxarcs,
        {{kk1+nn,kk2},{kk1,kk2+nn}}];
        ,{i,Length[subarcs[[center]]]}];
Do[{kk1,kk2}=subarcs[[center,i]];
      dist=Join[dist,
        {1-z[[kk1]]-z[[kk2]]-z[[center]],1-z[[kk1]]-z[[kk2]]-z[[center]]}];
        ,{i,Length[subarcs[[center]]]}];
Do[ result=shortestpathsym[subnodes[[center,i]],
                           subnodes[[center,i]]+nn,auxarcs,dist]; 
     If[! disconnected, result[[2]]=Mod[result[[2]]-1,nn]+1;
          result[[2]]=Drop[result[[2]],-1];
          If[result[[1]]<1-z[[center]],AppendTo[violations,result[[2]]]]
      ]
   ,{i,Length[subnodes[[center]]]}
];
If[violations!={},automaticaddwheelcuts[center]]
]

automaticaddwheelcuts[center_] := 
      Module[
{ordviol, sortedviol, marked, 
          essviol}, 
ordviol = 
          Map[goodorientation[#]  &, 
            violations]; 
sortedviol = 
          Union[Map[ord[#]  &, ordviol]]; 
marked = 
          Table[False, {Length[ sortedviol]}]; 
Do[
          If[Length[Union[sortedviol[[i]]]] <  
              Length[sortedviol[[i]]], 
            marked[[i]] = True], 
   {i, 
            Length[sortedviol]}]; 
essviol = Transpose[
              Select[Transpose[{sortedviol, 
                    marked}], !  #[[2]]  &]][[1]]; 
Do[nodeset = essviol[[i]]; 
          Print[ nodeset]; 
  If[
            Length[nodeset] == 
              3, 
       clique = 
              Union[nodeset, {center}]; 
      If[! MemberQ[fourcliques, 
                  clique], 
             AppendTo[fourcliques, 
                clique]; 
              addwheelcut[
                center]
           ]
     ,  
            addwheelcut[
              center]
    ];
        , {i, Length[essviol]}]]

addwheelcut[center_]:=Block[{a},Print["entering"];
 k=1;
col=Table[0,{m}];
Do[col[[nodeset[[i]]]]=1
   ,{i,Length[nodeset]}];
col[[center]]=(Length[nodeset]-1)/2;
If[z.col<=(Length[nodeset]-1)/2,
     Print["Though cut is ineffective the inequality was added"]];
Print[col];
n=n+1;
AT=Append[AT,col];
A=Transpose[AT];
c=Append[c,(Length[nodeset]-1)/2];
 ];

automaticaddcuts[v_] := 
      Module[
{ordviol, sortedviol, marked, 
          essviol}, 
ordviol = 
          Map[goodorientation[#]  &, 
            violations]; 
sortedviol = 
          Union[Map[ord[#]  &, ordviol]]; 
marked = 
          Table[False, {Length[ sortedviol]}]; 
Do[
          If[Length[Union[sortedviol[[i]]]] <  
              Length[sortedviol[[i]]], 
            marked[[i]] = True], 
   {i, 
            Length[sortedviol]}]; 
essviol = Transpose[
              Select[Transpose[{sortedviol, 
                    marked}], !  #[[2]]  &]][[1]]; 
Do[nodeset = essviol[[i]]; 
          Print[nodeset]; 
     addcircuitcut, {i, 
            Length[essviol]}]]

subset[a_, b_] := Union[a, b] == b

goodorientation[a_] := 
      Module[{p, b}, 
   b = 
          a; 
   p = Position[b, Min[a]][[1, 
              1]];  
    p0 = 
          Mod[p - 2, Length[a]] + 1; 
    p1 = 
          Mod[p, Length[a]] + 1; 
   If[
          a[[p1]] < a[[p0]], b, 
          Reverse[b]]
]

ord[a_] := 
      Module[{b, p}, 
   b = 
          Join[a, a];  
   p = Position[b, 
              Min[a]][[1, 1]];  
   Take[
          b, { p, p + Length[a] - 1}]]

preparezeta:=Module[{},
zfrac=Transpose[{Numerator[z],Denominator[z]}];
zstring=Map[Which[   #[[1]]==0
                   ," 0 "
                   , #[[2]]==1
                   ," 1 "
                   ,True
                   ,StringJoin[ToString[#[[1]]],"/",ToString[#[[2]]]]
                    ] & ,zfrac];
zfont=MapThread[FontForm,Transpose[Map[{#,{"cmbx10",12 }} &, zstring]]];
zlabels={zfont,nodes,displ};
zzz=Graphics[MapThread[Text,zlabels]];
]

findoddcircuit[z_]:=Module[{dist,result},
dist={};violations={};
Do[{kk1,kk2}=arcs[[i]];
      dist=Join[dist,
        {1-z[[kk1]]-z[[kk2]],1-z[[kk1]]-z[[kk2]]}];
        ,{i,Length[arcs]}];
Do[ result=shortestpathsym[i,i+nn,auxarcs,dist]; 
     If[! disconnected, result[[2]]=Mod[result[[2]]-1,nn]+1;
          result[[2]]=Drop[result[[2]],-1];
          If[result[[1]]<1,Print[result] ;AppendTo[violations,result[[2]]]]
      ]
   ,{i,nn}
]
]

shortestpathsym[s_,t_,arcs_,dist_]:=Block[
{rho,n,k,prev,mat,k1,k2,path,minrho,tree},
  n=Max[arcs]; 
  If[Length[arcs]!=Length[dist],
       Print["error in data for shortest path routine"];Abort[]];
  rho=Table[If[i==s,0,Infinity],{i,n}];
  disconnected=False;
  mat=Table[Infinity,{n},{n}];
  Do[k1=arcs[[i,1]];
     k2=arcs[[i,2]];
     mat[[k1,k2]]=dist[[i]];
     mat[[k2,k1]]=dist[[i]];
     ,{i,Length[arcs]}
     ];
  k=s;
  prev=Table[i,{i,n}];
  marked=Table[If[i==s,True,False],{i,n}];
  minrho=0;
  tree={};
  While[k!=t && minrho< Infinity,
      minrho=Infinity;
      Do[If[!marked[[i]]
            ,If[rho[[k]]+mat[[k,i]] < rho[[i]]
               ,rho[[i]]=rho[[k]]+mat[[k,i]];
                prev[[i]]=k;
               ];
             If[rho[[i]] < minrho
               ,minrho=rho[[i]];
                newk = i;
               ]
           ]
        ,{i,n}
        ];
      AppendTo[tree,{prev[[newk]],newk}];
      k=newk;
      marked[[k]]=True;
  ];
   If[minrho < Infinity
     ,path={t};
      k=t;
      While[k !=s,
          PrependTo[path,prev[[k]]];
          k=prev[[k]];
          ];
       {rho[[t]],path}
     ,disconnected=True
     ]
]

showedgecovering:=Block[{},
  smallbase={};
  Do[If[x[[i]]>0,AppendTo[smallbase,base[[i]]]],{i,Length[base]}];
  Show[pshow[Range[nn],2,"Bk"],show[arcs[[smallbase]],2,0,"R"],
        AspectRatio->Automatic,PlotRange->rangeplot];
]

show[a_,b_,c_,d_]:=Block[{b1,dd1,dd2,dd3},
           Which[b==1,b1=0.003,b==2,b1=0.005,
                 b==3,b1=0.008,b==4,b1=0.015];
           Which[d=="R",dd1=1;dd2=0;dd3=0;,
                 d=="B",dd1=0;dd2=0;dd3=1;,
                 d=="G",dd1=0;dd2=1;dd3=0;,
                 d=="Y",dd1=1;dd2=1;dd3=0.5;,
                 d=="Bk",dd1=0;dd2=0;dd3=0;
                 ];
           If[c==0,
          Graphics[{Thickness[b1],RGBColor[dd1,dd2,dd3],
          Table[Line[{nodes[[a[[i,1]]]],
                         nodes[[a[[i,2]]]]}
                       ]
                   ,{i,Length[a]}]}],
          Graphics[{Thickness[b1],Dashing[{0.01}],RGBColor[dd1,dd2,dd3],
          Table[Line[{nodes[[a[[i,1]]]],
                         nodes[[a[[i,2]]]]}
                       ]
                   ,{i,Length[a]}]}]
                   ]]

pshow[a_,b_,d_]:=Block[{b1,dd1,dd2,dd3},
           Which[b==1,b1=0.01,b==2,b1=0.022,
                 b==3,b1=0.035,b==4,b1=0.05,
                 b==0,Null];
           Which[d=="R",dd1=1;dd2=0;dd3=0;,
                 d=="B",dd1=0;dd2=0;dd3=1;,
                 d=="G",dd1=0;dd2=1;dd3=0;,
                 d=="Y",dd1=1;dd2=1;dd3=0.5;,
                 d=="Bk",dd1=0;dd2=0;dd3=0;
                 ];
           If[b==0,
                Graphics[{ RGBColor[dd1,dd2,dd3],
                            Table[{PointSize[0.01+z[[a[[i]]]] 0.05],
                                Point[nodes[[a[[i]]]]]},
                                {i,Length[a]}]
                                }
                           ]
                   ,
            Graphics[{PointSize[b1],RGBColor[dd1,dd2,dd3],
          Table[Point[nodes[[a[[i]]]]]
                   ,{i,Length[a]}]}]
              ]  ]

checkhole[a_]:=Block[{fail,arc},
      fail=False;
      {firstcheck,arc}=checkarcs[a];
      If[firstcheck,
         If[OddQ[Length[a]]
            ,{holetest,arc}=yeshole[a];
              fail=!holetest;
              If[fail,Print["Not a hole ! ",arc," is a chord"]];
             ,Print["Even hole !"];
              fail=True
            ]
            ,
            Print["Arc ",arc," does not exist !"];
            fail=True
          ];
      !fail
]

checkarcs[a_]:=Block[{b,k,test},
     b=Append[a,a[[1]]];
     k=1;
     test=True;
     While[k<=Length[a] && test,
           test=MemberQ[arcs,{b[[k]],b[[k+1]]}] ||
                MemberQ[arcs,{b[[k+1]],b[[k]]}];
           k++
          ];
      {test,{b[[k-1]],b[[k]]}}
]

yeshole[a_]:=Block[{k1,k2,test},
     k1=1;k2=3;
     test=False;
     While[!test && (k2<=Length[a]-1)
                ,
                 test=MemberQ[arcs,{a[[k1]],a[[k2]]}] ||
                      MemberQ[arcs,{a[[k2]],a[[k1]]}];
                  k2++
                  ];
            k1++;
     While[!test && (k1<Length[a]-1)
          ,k2=k1+2;
           While[!test && (k2<=Length[a])
                ,
                 test=MemberQ[arcs,{a[[k1]],a[[k2]]}] ||
                      MemberQ[arcs,{a[[k2]],a[[k1]]}];
                  k2++
                  ];
            k1++;
          ];
      {!test,{a[[k1-1]],a[[k2-1]]}}
]

addholecut:=Block[{a},
If[checkhole[nodeset],
k=1;
col=Table[0,{m}];
Do[col[[nodeset[[i]]]]=1
   ,{i,Length[nodeset]}];
If[z.col<=(Length[nodeset]-1)/2,
     Print["Though cut is ineffective the inequality was added"]];
n=n+1;
AT=Append[AT,col];
A=Transpose[AT];
c=Append[c,(Length[nodeset]-1)/2];
,
Print["Cannot add inequality !"];
]

];

addcircuitcut:=Block[{a},
 k=1;
col=Table[0,{m}];
Do[col[[nodeset[[i]]]]=1
   ,{i,Length[nodeset]}];
If[z.col<=(Length[nodeset]-1)/2,
     Print["Though cut is ineffective the inequality was added"]];
n=n+1;
AT=Append[AT,col];
A=Transpose[AT];
c=Append[c,(Length[nodeset]-1)/2];
 ];

branchto[node_,v_]:=Block[{a},
 k=1;
col=Table[0,{m}];
If[v==0,col[[node]]=1,col[[node]]=-1];
n=n+1;
AT=Append[AT,col];
A=Transpose[AT];
c=Append[c,-v];
 ];

checkclique[a_]:=Block[{k1,k2,test},
      test=True;
      k1=1;
      While[test && (k1<Length[a]),
            k2=k1+1;
            While[test  &&  (k2<=Length[a]),
                 test=MemberQ[arcs,{a[[k1]],a[[k2]]}] ||
                       MemberQ[arcs,{a[[k2]],a[[k1]]}];
                  k2++;
                   ];
            k1++;
           ];
       {test,{a[[k1-1]],a[[k2-1]]}}
   ]

addcliquecut:=Block[{a,firstcheck,arc},
{firstcheck,arc}=checkclique[nodeset];
If[firstcheck,
k=1;
col=Table[0,{m}];
Do[col[[nodeset[[i]]]]=1
   ,{i,Length[nodeset]}];
If[z.col<=1,Print["Though cut is ineffective the inequality was added"]];
n=n+1;
AT=Append[AT,col];
A=Transpose[AT];
c=Append[c,1];
,
Print["Not a clique ! ","arc ",arc," is missing"]
]
];

addcut[a_, b_]:=Block[{  },  
 
        n = n + 1; 
AT = 
          Append[AT, a]; 
A = 
          Transpose[AT]; 
c = 
          Append[c, b]]

computex:=Block[{},
optimal=False;
infeasible=False;
While[! optimal && ! infeasible,
      If[newbase,
         B = Transpose[Transpose[A][[ base ]]];
         W=Inverse[B];
         x = W. b;
         cb = c[[ base ]];Print[N[cb.x]];
         z = cb.W;
         newbase=False];
      cr = c - z . A;
      Block[{i},
            For[i=1,
                (i<=n) && (cr[[i]]>=-0.001),
                i++,Null];
             If[i==n+1,optimal=True,p=i]   
           ] ;
      If[optimal,Null ,
         ap=Transpose[A][[p]];
         pivot;  ];
      ];
 ];

pivot:=Block[{i,ratio,min},
      y=W.ap;
      min = 1000000;
      For[i=1, i<=m, i++,
          If[ y[[i]]>0,
              ratio = x[[i]]/y[[i]];
              If[ratio < min,
                 min = ratio; q = i]
            ]  
      ];
      base[[q]]=p;
      newbase=True;
      infeasible=(min==1000000);
];

initialize:=Block[{AT},
     AT=Table[0,{i,nn},{j,na}];
     k=1;
     Do[AT[[i,k]]=1;AT[[j,k]]=1;k++;
        ,{i,nn-1},{j,i+1,nn}];
     c=Join[Table[1,{nn}],Table[-1,{nn}],Table[0,{na}]];
     A=Transpose[Join[AT,-AT,-IdentityMatrix[na]]];
     base=Table[2 nn+i,{i,na}];
     b=-bb;
     newbase=True;
     n=2nn+na;
     m=na;
     i=.;arclist={};
     k=1;
     Do[AppendTo[arclist,{i,j}];
      k++,
        {i,nn-1},{j,i+1,nn}];
      ];

Data

User provided data

nodes is the list of node coordinates
arcs is the list of arcs

Set #1

2 clique cuts

nodes={{64, 20}, {36, 47}, {50, 60}, {90, 45},
 {10, 30}, {20, 0}, {70, 44}, 
   {43, 22}, {80, 20}, {40, 0}};

nn = Length[nodes]

labeldispl = Table[{-2, -2}, {nn}];

arcs={{1, 2}, {1, 3}, {1, 7}, {1, 8}, {1, 10}, {2, 3}, {2, 5},
 {2, 6},  {2, 7}, {2, 8}, {3, 4}, {3, 7},
  {3, 8}, {4, 7}, {4, 9}, {6, 10}, {7, 9}};

Set #2

(intermediate instance)
hint: use 
cliques  {2,6,11}; {9,4,11}; {7,4,12}
oddhole  {1,8,12,7,9,5,3}

nodes={{80, 98}, {80, 44}, {20, 50}, {17, 63}, {59, 4},
     {63, 30}, {27, 92}, 
   {51, 90}, {82, 33}, {9, 22}, {50, 52}, {78, 70}};

nn = Length[nodes]

labeldispl = Table[{-2, -2}, {nn}];

arcs={{1, 3}, {1, 6}, {1, 8}, {1, 10}, {1, 11}, {2, 6}, {2, 8}, 
    {2, 10}, 
   {2, 11}, {3, 5}, {3, 6}, {4, 7}, {4, 9}, {4, 11}, {4, 12}, {5, 9}, 
   {5, 10}, {6, 11}, {6, 12}, {7, 9}, {7, 12}, {8, 12}, {9, 11}};

Set #3

trivial, only 3 clique cuts

nodes={{0,0},{1,0},{2,0},{3,0},
          {0.5,1},{1.5,1},{2.5,1},
             {1,2},{2,2},
                 {1.5,3}};

nn = Length[nodes]

labeldispl = Table[{-2, -2}, {nn}];

arcs={{1,2},{2,3},{3,4},{5,6},{6,7},{8,9},
      {1,5},{5,8},{8,10},{2,6},{6,9},{3,7},
      {10,9},{9,7},{7,4},{8,6},{6,3},{2,5}};

Set #4

nodes={{0,0}, {1,0}, {2,0},       {4,0},      {6,0},
              {1,1},
       {0,2},        {2,2}, {3,2},       {5,2},{6,2},
       {0,3}, {1,3}, {2,3},       {4,3},      {6,3},
                     {2,4},
       {0,5},        {2,5},                   {6,5}};

nn = Length[nodes]

labeldispl = 
        Table[{-2, -2}, {nn}];
labeldispl[[15]] = {0, -2}; 
    labeldispl[[6]] = {0, -2};
labeldispl[[7]] = {-2, -3}; 
    labeldispl[[10]] = {0, -2};
labeldispl[[1]] = {-2, -3}; 
    labeldispl[[4]] = {3, -2}; 
    labeldispl[[3]] = {-3, -2};

arcs={{1,2},{2,3},{3,4},{4,5},{1,6},{2,6},{3,6},
   {1,7},{6,7},{3,8},{3,9},{4,9},{4,15},{4,10},{4,11},
   {5,11},{6,8},{7,8},{7,13},{7,12},{8,13},{8,14},{9,14},
   {10,15},{11,15},{11,16},{12,13},{12,18},{13,14},{13,18},
   {13,19},{14,15},{14,17},{15,17},{15,19},{15,20},{15,16},
   {16,20},{17,19},{18,19},{19,20},{10,11}};

Set #5

nn=7;

nodes=Table[5 {Sin[2 N[Pi]/nn i],Cos[2 N[Pi]/nn i]},{i,0,nn-1}];

arcs={{1,2},{2,3},{3,4},{4,5},{5,6},{6,7},{1,7},
     {1,3},{3,5},{5,7},{2,7},{2,4},{4,6},{1,6}};

Set #6

nodes={{0,0}, {1,0}, {2,0},       {4,0}, 
       {0,1},        {2,1},
              {1,2},        {3,2},{4,2},  
       {0,3},        {2,3},       {4,3}};

nn = Length[nodes]

labeldispl = Table[{-2, -2}, {nn}];

arcs={ {1,2},{1,5},{2,3},{2,5},{3,4},{3,6},{4,6},{4,9},
{5,6},{5,7},{5,10},{6,7},{6,8},{6,9},{6,11},
{7,10},{7,11},{8,9},{8,11},{9,11},{9,12},{10,11},{11,12}};

Set #7

nodes={{0,0}, {0,1}, {0,2},            
              {1,1},
       {2,1},        {2,2},          
       {3,0}, {3,1}, {3,2},               
                     {4,2},
       {5,0},        {5,2}                   };

nn = Length[nodes];

labeldispl = 
        Table[{-2, -2}, {nn}];
labeldispl[[6]] = {0, -2};
labeldispl[[7]] = {-2, -3}; 
    labeldispl[[10]] = {0, -2};
labeldispl[[1]] = {-2, -3}; 
    labeldispl[[4]] = {3, -2}; 
    labeldispl[[3]] = {-3, -2};

arcs={{1, 2}, {2, 3}, {1, 4}, {2, 4}, {3, 4}, 
{1, 5}, {4, 5}, {3, 6}, {4, 6}, {1,7},
 {5, 6}, {5, 8}, {5, 7}, {6, 8}, 
 {6, 9}, {7, 8}, {7, 11}, {8, 9}, {8, 11}, 
 {8, 12}, {9, 10}, {10, 12}, {11, 12}}

Data preprocessing (must be done for any data set)

nn=Length[nodes];
na=Length[arcs];

auxarcs={};
Do[{kk1,kk2}=arcs[[i]];
     auxarcs=Join[auxarcs,
        {{kk1+nn,kk2},{kk1,kk2+nn}}];
        ,{i,Length[arcs]}];

adjtable = Table[{}, {nn}];
Do[  {k1, k2} = 
          arcs[[i]]; 
         AppendTo[
          adjtable[[k1]], 
          k2]; 
         AppendTo[
          adjtable[[k2]], k1];, {i, 
          Length[arcs]}];
subarcs = {}; subnodes = {};
Do[ 
        AppendTo[
          subarcs, 
           Select[arcs,  
            MemberQ[adjtable[[i]], #[[1]]]   &&  
                                     MemberQ[
                  adjtable[[i]], #[[2]]]   & 
                        ]
                    ]; 
       AppendTo[
          subnodes, 
          Union[Flatten[
              subarcs[[-1]]]]]
        , {i, nn}];

minx=Min[Transpose[nodes][[1]]];
miny=Min[Transpose[nodes][[2]]];
maxx=Max[Transpose[nodes][[1]]];
maxy=Max[Transpose[nodes][[2]]];

rx=maxx-minx;
ry=maxy-miny;

rangeplot={{minx-.2 rx,maxx+.2 rx},{miny-.2 ry,maxy+.2 ry}}

nodelabelfont = 
        MapThread[
          FontForm, 
Transpose[
            Map[{#, {"Avant Garde", 12}} &, 
              Range[nn]]]];
nodelabels = {nodelabelfont, nodes, 
          labeldispl};
aaa = 
        Graphics[{PointSize[0.02], 
            Map[Point, nodes], 
MapThread[Text, 
              nodelabels ]}];

(*aaa=Graphics[Table[{PointSize[0.0],Point[nodes[[i]]],
      Text[ToString[i],nodes[[i]],{-2,-2}]},{i,Length[nodes]}]];*)

displ=Table[{2,2},{i,Length[nodes]}];

displ[[6]] = {-3, 2}; displ[[9]] = {-2, 2}; 
    displ[[10]] = {0, 2};
displ[[12]] = {-2, 0}; displ[[3]] = {2, 0}; 
    displ[[1]] = {2, 0}; 
    displ[[2]] = {2, 0};
displ[[4]] = {-1, 2.5}; displ[[5]] = {0, 2.5};  
    displ[[8]] = {-3, 0};
displ[[11]] = {-2, 0}; 
    displ[[7]] = {2.5, -1};

ccc=pshow[Range[nn],2,"Bk"];

bbb=show[arcs,1,0,"B"];

labeldispl = 
        Table[{-2, -2}, {nn}];
labeldispl[[6]] = {0, -2}; labeldispl[[9]] = {0, -2}; 
    labeldispl[[12]] = {0, -2};
labeldispl[[7]] = {-2, -3}; 
    labeldispl[[10]] = {0, -2}; 
    labeldispl[[11]] = {1, -2};
labeldispl[[1]] = {-2, -3}; 
    labeldispl[[4]] = {0, -2}; 
    labeldispl[[3]] = {0, -2};

Show[aaa,bbb,ccc,AspectRatio->Automatic,PlotRange->rangeplot];

Randomly generated data

Graphs are generated randomly; for each undirected pair (i,j) 
the arc (i,j) is generated with probability d/(nn-1), with
nn the  number of nodes. The user should provide nn,d and seed.

Don't do any computation if the graph is disconnected !

interesting values: 
   seed=60 , nn=10,  d=3;   (very quick)
   seed=30 , nn=12,  d=3;   (intermediate)
   seed=20 , nn=24,  d=5;   (tough)

nn=12;
seed=20;
d=5;

SeedRandom[seed];
prob=d/(nn-1);
arcs={}; 
Do[If[Random[]<prob,
      AppendTo[arcs,{i,j}];
      ],
      {i,1,nn-1},{j,i+1,nn}];

nonsinglenodes = Union[Flatten[arcs]];

isolatednodes = Complement[Range[nn], nonsinglenodes];

Do[
        AppendTo[
          arcs, {isolatednodes[[i]], 
nonsinglenodes[[Random[
                  Integer, {1, 
                    Length[nonsinglenodes]}]]]}], 
{i, 
          Length[isolatednodes]}];
na = Length[arcs];

auxarcs={};
Do[{kk1,kk2}=arcs[[i]];
     auxarcs=Join[auxarcs,
        {{kk1+nn,kk2},{kk1,kk2+nn}}];
        ,{i,Length[arcs]}];

adjtable = Table[{}, {nn}];
Do[  {k1, k2} = 
          arcs[[i]]; 
         AppendTo[
          adjtable[[k1]], 
          k2]; 
         AppendTo[
          adjtable[[k2]], k1];, {i, 
          Length[arcs]}];
subarcs = {}; subnodes = {};
Do[ 
        AppendTo[
          subarcs, 
           Select[arcs,  
            MemberQ[adjtable[[i]], #[[1]]]   &&  
                                     MemberQ[
                  adjtable[[i]], #[[2]]]   & 
                        ]
                    ]; 
       AppendTo[
          subnodes, 
          Union[Flatten[
              subarcs[[-1]]]]]
        , {i, nn}];

arcs;

complarcs = 
        Complement[Flatten[Table[{i, j}, {i, nn - 1}, {j, i + 1, nn}], 1], 
          arcs];

nodes=Table[5 {Sin[2 N[Pi]/nn i],Cos[2 N[Pi]/nn i]},{i,0,nn-1}];

rangeplot=5.4{{-1.2,1.2},{-1.2,1.2}};

aaa=Graphics[Table[{PointSize[0.0],Point[nodes[[i]]],
      Text[ToString[i],nodes[[i]],0.6{-nodes[[i,1]],-nodes[[i,2]]}]},{i,Length[nodes]}]];

displ=Table[0.8{-nodes[[i,1]],-nodes[[i,2]]},{i,Length[nodes]}];

ccc=pshow[Range[nn],2,"Bk"];

bbb=show[arcs,1,0,"B"];

zzz = show[complarcs, 1, 0, "R"]

Show[aaa,bbb,ccc,AspectRatio->Automatic,PlotRange->rangeplot];

Show[aaa,zzz,ccc,AspectRatio->Automatic,PlotRange->rangeplot];

Phase I

This finds an initial feasible solution to the dual problem.
x is the dual variable associated to the arcs. The artificial objective is displayed.  Phase I is completed when the objective goes to 0.
At this point x (if integral) is an edge covering such that its cardinality plus the number of nodes with degree greater than 1 is not greater than the number of nodes

n=na+2 nn;m=nn;
AT={};
Do[k1=arcs[[j,1]];k2=arcs[[j,2]];
   AppendTo[AT,Table[If[i==k1 || i==k2,1,0],{i,nn}]]
   ,{j,na}];
AT=Join[AT,-IdentityMatrix[nn],IdentityMatrix[nn]];
A=Transpose[AT];
c=Table[If[j>na+nn,1,0],{j,n}];
b=Table[1,{nn}];
base=Table[j,{j,na+nn+1,n}];
newbase=True;

computex

x

AT=Drop[AT,-nn];
A=Transpose[AT];
n=na+nn;
c=Table[If[j<=na,1,0],{j,n}];
newbase=True;

(*showedgecovering*)

Phase II

Reoptimization

LP solution

Solution of the dual of the fractional node packing
system plus added cut inequalities

computex

z

Display

Want to see node values? (write True or False, possibly adjust displ if True)

seenodevalues=False;

regset={};
fractset={};
zeroset={};
Do[Which[z[[i]]==1,
         AppendTo[regset,i],
         z[[i]]==0,
         AppendTo[zeroset,i],
         True,
         AppendTo[fractset,i]
         ]
   ,
  {i,Length[base]}];
  Print["integral nodes  = ",regset];
  Print["fractional nodes= ",fractset];

ccc=pshow[fractset,0,"R"];

ddd=pshow[regset,4,"G"];

eee=pshow[zeroset,2,"B"];

If[seenodevalues==True,
   preparezeta;
   Show[aaa,bbb,ccc,ddd,eee,zzz,AspectRatio->Automatic,
       PlotRange->1.1 rangeplot];
       ,
   Show[aaa,bbb,ccc,ddd,eee, AspectRatio->Automatic,
     PlotRange->1.03 rangeplot];
 ]

Cutting

Semiautomatic

getting odd circuits

findoddcircuit[z]

automatic adding odd circuit cuts (must be preceeded by findoddcircuit[z]

automaticaddcuts[violations]

getting odd wheels (to be used only if there are no violated odd hole inequalities)

fourcliques = {};
Do[If[subarcs[[i]] ≠ {}, Print["center ", i]; 
        findoddwheel[i, z]], {i, nn}]

Manual

Clique inequalities

nodeset={1,2,6};
addcliquecut;

Odd hole inequalities

nodeset={7,8,14,17,19,18,12};
addholecut

Finding antihole inequalitites

specnodes = 
      Select[Range[nn],  
        z[[#]] > 0  &&  z[[#]] < 1  &]

complarcs = {};
Do[ 
      ii = specnodes[[i]]; 
      jj = 
        specnodes[[j]]; 
      If[! MemberQ[
            arcs, {ii, jj}],  
        AppendTo[complarcs, {ii, jj}]]  , 
{i, 
        Length[specnodes] - 1}, 
 {j, i + 1, 
        Length[specnodes]}]

complarcs

fff = show[complarcs, 1, 0, "B"];

This is the complement of the subgraph restircted to the fractional nodes

Show[aaa, fff, ccc, ddd, eee, AspectRatio -> Automatic, 
        PlotRange -> 1.03 rangeplot];

Identify  a odd hole with at least five nodes in the graph

antihole = {4, 15, 9, 7, 13, 10}

checking for chords

oddholearcs =  
      Map[Sort,   Transpose[ {antihole, 
            RotateRight[antihole]}]]

Do[If[MemberQ[complarcs, oddholearcs[[i]]], Nil, 
        Print[oddholearcs[[i]], " not a complementary arc"]], {i, 
        Length[antihole]}]

chords = 
      Complement[
        Select[complarcs, 
          MemberQ[antihole,  #[[1]]]  && 
                                                              MemberQ[antihole,  #[[2]]]   &], 
        oddholearcs]

If OK go on

antiholearcs = 
      Select[arcs, 
        MemberQ[antihole,  #[[1]]]  && 
                                                              MemberQ[
              antihole,  #[[2]]]   &]

Apply[Plus, Map[z[[#]]  &,  antihole]]

The odd hole inequality is violated  if the above number is greater than 2

a = Table[0, {nn}]

Do[a[[antihole[[i]]]] = 1, {i, Length[antihole]}]

a