model columngen.mod;
data columngen.dat;

      
for {k in K} {   # compute the cost of an edgeset
  let Nodes[k] := {};
  for {(i,j) in Edgeset[k]}  {
     let Nodes[k] := Nodes[k] union {i} union {j};
#     let a[i,j,k] := 1;
  }
}
let {k in K} c[k] := card(Nodes[k]);
let {(i,j) in E, k in K} a[i,j,k] := card({(i,j)} intersect Edgeset[k]);


for {counter in 1 .. 100} {

  problem Master;

  printf "Solving restricted LP master problem with the following rings:\n";
  display Edgeset;

  option relax_integrality 1;
  option solver cplex;
  solve;
  for {k in K: z[k] > 0} {
     printf "z[%d] = %f\n",k,z[k];
     display Edgeset[k];
  }

 printf "Dual variables:\n";
 display edgecover;


  printf "Solving the pricing problem:\n";
  let {(i,j) in E} pi[i,j] := edgecover[i,j];
  problem Pricing;

  expand reduced_cost;
  expand maxlength;
  option relax_integrality 0;
  option solver cplex;
  solve;
 # Stop with the column generation if the optimal
 # solution of the pricing problem is close to 0

  if reduced_cost >= -1e-10 then
    break;

  display x;
  display y;
  printf "New edgeset\n";
  display {(i,j) in E: y[i,j] == 1};

  let h := card(K) + 1;
  let K := K union {h};
  let Edgeset[h] := {(i,j) in E: y[i,j] == 1};
  let Nodes[h] := {};
  for {(i,j) in Edgeset[h]} {
     let Nodes[h] := Nodes[h] union {i} union {j};
#     let a[i,j,h] :=1;
  }
  let c[h] := card( Nodes[h] );

  let {(i,j) in E} a[i,j,h] := card({(i,j)} intersect Edgeset[h]);
} # end for {counter in 1..100}


printf "Solve the IP:\n";
display Edgeset;
problem Master;
option relax_integrality 0;
expand;
solve;

for {k in K: z[k] > 0} {
 printf "z[%d] = %f\n",k,z[k];
 display Edgeset[k];
}