/************************************************************************************

 * The model shows a Benders Decomposition for a location-transportation problem.   *

 * The original MIP is decomposed into two problems.                                *

 * The subproblem is using the dual formulation.            			    *

 * Version: 1.0, April 10th, 2005, Tested with CPLEX v9.0			    *					

 * For consistency, this is built using an equivalent AMPL model cfl.benders.mod    *

 ************************************************************************************/










#include <ilcplex/ilocplex.h>



ILOSTLBEGIN



typedef IloArray<IloNumArray> TwoDMatrix;






int

main(int argc, char **argv) {

 IloEnv env;



 try {



	 IloInt i=0, j=0, k=0;

	 IloInt ORIG=0, DEST=0;

	 IloInt MaxCut=0;               //Maximum number of cuts 


	 IloNumArray supply(env), demand(env);

	 IloNumArray fix_cost(env);
	 TwoDMatrix var_cost(env);







///////////////// DATA FILE READING ////////////////////////////////

	  const char* filename  = "cfl_benders.dat";

          if (argc > 1)

             filename = argv[1];

          ifstream file(filename);

          if (!file) {

             cerr << "ERROR: could not open file '" << filename

                  << "' for reading" << endl;

             cerr << "usage:   " << argv[0] << " <file>" << endl;

             throw(-1);

          }



       	  file >> MaxCut >> supply >> demand >> fix_cost >> var_cost ;

	  ORIG = supply.getSize();
	  DEST = demand.getSize();




 	  IloBool consistentData = (fix_cost.getSize() == ORIG && var_cost.getSize() == ORIG);

    	  if (!consistentData) {

 	      cerr << "ERROR: data file '"

                   << filename << "' contains inconsistent data" << endl;

              throw(-1);

	   }







/////////////////// DECISION VARIABLES WITH NAMES  /////////////////////////////


          ////////DECISION VARIABLES AND PARAMETERS FOR SUBPROBLEM///////////////

 	   IloNumArray build(env, ORIG);      // = 1 iff facility built at i



 	   IloNumVarArray Supply_Price(env, ORIG, -IloInfinity, 0, ILOFLOAT);

	   IloNumVarArray Demand_Price(env, DEST, -IloInfinity, IloInfinity, ILOFLOAT);




	  ////////DECISION VARIABLES AND PARAMETERS FOR MASTER PROBLEM///////////

	   IloInt nCut = 0;

	   IloNumArray cut_type(env, MaxCut);  //= 0 iff subproblem has optimal solution; 
		                               //= 1 iff subproblem is unbounded



	   TwoDMatrix supply_price(env, MaxCut);
	   TwoDMatrix demand_price(env, MaxCut);

	   for(i=0; i<MaxCut; i++){

		supply_price[i] = IloNumArray(env, ORIG);
		demand_price[i] = IloNumArray(env, DEST);

	   }



  	   IloBoolVarArray  Build(env, ORIG);
	   IloNumVar Max_Ship_Cost(env, 0, IloInfinity, ILOFLOAT);




////////////DEVELOP GENERIC MODEL //////////////////////////


	    //////SET MASTER PROBLEM///////////////////////////
		
	    IloModel model_master(env);

  	    IloExpr Objective_master(env);

	    Objective_master += Max_Ship_Cost;

	    for(i=0; i<ORIG; i++){

		Objective_master += fix_cost[i] * Build[i];

	    }

            model_master.add(IloMinimize(env, Objective_master));


	    IloCplex cplex_master(env);



            ///////SET SUBPROBLEM (DUAL FORMULATION)//////////////////////

	    IloModel model_sub(env);



	    IloObjective Objective_sub(env);

 	    model_sub.add(Objective_sub);




 	    for(i=0; i<ORIG; i++){

		for(j=0; j<DEST; j++){

			model_sub.add(Supply_Price[i] + Demand_Price[j] <= var_cost[i][j]);

		}

	     }


	     IloCplex cplex_sub(model_sub);

	     IloNum eps = cplex_sub.getParam(IloCplex::EpInt);





//////////BEGIN ITERATIONS/////////////////////////////////


		IloNum GAP = IloInfinity;
		IloNum max_ship_cost = 0;
		IloNum Dual_Ship_Cost =0;


		for(i=0; i<ORIG; i++){

			build[i] = 1;
		}

		while( nCut < MaxCut ){

			env.out() << " ITERATION "<< nCut+1 << endl << endl; 


			IloExpr sub_obj(env);

			for(i=0; i<ORIG; i++){
				
				sub_obj += -supply[i] * build[i] * Supply_Price[i];

			}
			
			for(j=0; j<DEST; j++){

				sub_obj += -demand[j] * Demand_Price[j];

			}


			Objective_sub.setExpr(IloMinimize(env, sub_obj));

			

			cplex_sub.setParam(cplex_sub.PreInd, 0);   //Disable presolve, otherwise, if dual is infeasible, 
                                                                   //we don't know if prime is unbounded or infeasible


			cplex_sub.solve();

			
			if( cplex_sub.getCplexStatus() == CPX_STAT_UNBOUNDED){ //subproblem is unbounded

				cut_type[nCut] = 1;


				cplex_sub.getRay( supply_price[nCut], Supply_Price);
				cplex_sub.getRay( demand_price[nCut], Demand_Price);

			
			}
			else{

				
				
				Dual_Ship_Cost = - cplex_sub.getObjValue();

				
				env.out()<<endl<< "Dual_Ship_Cost = "<<Dual_Ship_Cost<<endl;

				
				if( Dual_Ship_Cost <= max_ship_cost + eps){
					
					break;
				}

				if( GAP + eps > Dual_Ship_Cost - max_ship_cost){

					GAP = Dual_Ship_Cost - max_ship_cost;

				}

				env.out()<< "GAP = "<< GAP << endl;

				cut_type[nCut]=0;

				cplex_sub.getValues( supply_price[nCut], Supply_Price);
				cplex_sub.getValues( demand_price[nCut], Demand_Price);

			}


			// Add constraints to master problem

        		IloExpr lhs(env);

			for(i=0; i<ORIG; i++){

				lhs += supply_price[nCut][i] * supply[i] * Build[i];

			}

			IloNum rhs = 0;

			for(j =0; j<DEST; j++){

        			rhs += demand_price[nCut][j] * demand[j];

			}

			if( cut_type[nCut] == 1){

				model_master.add( lhs <= -rhs);

			}
			else{

				model_master.add( Max_Ship_Cost >= lhs + rhs); 

			}
			
			nCut++;

			env.out()<< endl << " RE-SOLVING MASTER PROBLEM "<<endl <<endl;



			cplex_master.extract(model_master);

			if( !cplex_master.solve () ){
				
				env.error() << "Failed" << endl;
				throw(-1);
			
			}  

			


			cplex_master.getValues( build, Build);
			max_ship_cost = cplex_master.getValue(Max_Ship_Cost);

			env.out()<< endl<< "Max_Ship_Cost = " << max_ship_cost<<endl;
			env.out()<<endl<<"Build [*] = " << endl << build<<endl;

		}

                if(nCut == MaxCut){
		     
		  env.out()<<endl<<"Reach limit of Cuts."<<endl;
                  
                }
                else{

		  env.out()<<endl<<endl<<"Optimal Value (Master Problem) = "<<cplex_master.getObjValue()<<endl;
                  env.out()<<endl<<"Optimal Value of Subproblem (Dual_Ship_Cost): "<<Dual_Ship_Cost<<endl<<endl;
                  
                  for(i=0; i<ORIG; i++){

		    if(build[i] >1-eps){

		      env.out()<<endl<<"Build["<<i+1<<"] = 1 ;"<<endl;

                    }
		  }
                }

		model_master.end();
		model_sub.end();
		cplex_master.end();
		cplex_sub.end();









 }

 catch(IloException &e) {

	env.out() << "ERROR: " << e << endl;

 }

 catch(...){

	env.out() << "Unknown exception" << endl;

 }



 env.end();



return 0;

}