/************************************************************
*
* MST solver using integer programming
*
*************************************************************/

// Prepare an open array to hold the components generated during the
// algorithm. Each of these define one or more cuts.


Open setof(int) componentset[1..0];

Model m("mst.mod", "gr17.dat");
//Model m("mst.mod", "ch150.dat");
//Model m("mst.mod","gr6.dat");
//Model m("mst.mod", "data\gr21.dat");
//Model m("mst.mod", "data\gr24.dat");
//Model m("mst.mod", "data\ulysses16.dat");
//Model m("mst.mod", "data\ulysses22.dat");
//Model m("mst.mod", "data\berlin52.dat");
//Model m("mst.mod", "data\ch130.dat");
//Model m("mst.mod", "data\ch150.dat");
//Model m("mst.mod", "data\dantzig42.dat");
//Model m("mst.mod", "data\fri26.dat");
//Model m("mst.mod", "data\gr48.dat");
//Model m("mst.mod", "data\hk48.dat");
//Model m("mst.mod", "data\pr136.dat");
//Model m("mst.mod", "data\pr76.dat");
//Model m("mst.mod", "data\rat99.dat");

//Data and statistics//
int n               := m.n;  // # vertices in G
range vertices      1..n;
int       component[1..n];
int       newcomponents;
int new;
int z;
int mark[1..n];
int scanned[1..n];
int iter :=0;
float time := 0.0;




//Main loop
repeat {

       // Solve MIP and retrieve results
       m.reset();
       m.solve();
       cout << "Objective =" << m.objectiveValue()<< endl;
       iter := iter+1;
       time := time+m.getTime();


       // Initialize connected component data
       newcomponents := 0;
       setof(int) firstcomp :={};
       setof(int) covered :={};

       // find all connected components

       // initialize all vertices  as unmarked and unscanned

       forall (i in vertices) mark[i] := n;
       forall (i in vertices) scanned[i] := 0;
       int nbscanned := 0;  // no vertices scanned

       repeat {

         //pick smallest indexed vertex not scanned yet
         int q:= min (i in vertices : scanned[i]=0) i;
         
          //Find component containing the vertex q
         
         
         new := 1;
         int p := 0;
         mark[q] := 0;
         repeat {
           p := p+1;
           z := min (i in vertices : scanned[i]=0) mark[i];
           q := min (i in vertices : (scanned [i]=0  & mark[i]=z)) i;
           
           
           component[p] := q;
           scanned[component[p]] :=1; //mark vertex as scanned
           
           nbscanned:= nbscanned+1;
           
           forall (i in vertices) forall (j in vertices : m.y[i,j]=1){
             if (mark[j] =n & component[p]=i) then{
                   mark[j]:=mark[component[p]]+1;
                   new :=new+1; 
                   
                 } 
             else 
                 if (mark[i]=n & component[p]=j) then {
                   mark[i]:= mark[component[p]]+1;
                   new:=new+1;
                 }
           }       
             
           
         } until (nbscanned = new \/ nbscanned = n);
         
          setof(int) comp := {component[i] | i in 1..p};
          
          covered := covered union comp;
          cout << "no elements covered =" << card(covered)<<endl;
         
          // Add component

          componentset.addh();
          newcomponents:=newcomponents+1;
          componentset[componentset.up]:=comp;
       } until card(covered) = n;          
              
       cout << "objective =" << m.objectiveValue() <<",";
       cout << newcomponents << " new component constraint(s), ";
       cout << m.getTime() << " sec " << endl;
         

} until newcomponents = 1;  //found a spanning tree

// print results

cout << endl << "*****************************************" << endl;
cout << iter << " iterations" << endl;
cout << componentset.size << " components resolved" << endl;
cout << "Optimal tree: ";

forall (i in vertices) forall (j in vertices  : (i< j & m.y[i,j]=1))
  cout <<  " < " << i << "," << j << ">  ";
cout << endl;
cout << "Total computation time: " << time << " sec" << endl;