// Dec. 2, 2012
// Coded by Minjiao Zhang, based on Dr. Kucukyavuz and Saumya's work
// Dynamic probablisitic Lot-Sizing with fixed cost

#include <ilcplex/ilocplex.h>
#include <stdlib.h>
#include <math.h>
#include <cstdlib> 
#include <time.h>

#define eps 0.00001
#define eps1 0.00001

ILOSTLBEGIN

IloInt n, m;
IloNum tau;
IloInt nbMix=0;
IloInt nbCmix=0;

IloInt nbIteration=0;

typedef IloArray<IloNumArray> NumMatrix;
typedef IloArray<NumMatrix> NumMatrix3d;
typedef IloArray<IloNumVarArray> NumVarMatrix;
typedef IloArray<IloBoolVarArray> BoolVarMatrix;

typedef IloArray<IloExpr> ExprArray;
typedef IloArray<ExprArray> ExprMatrix;

IloNum comp(IloNum f1, IloNum f2)	// for function f_1^2
{
	if(f1>f2-eps1)
	{
		return f1-f2;
	}
	else
	{
		return f1-f2+1;
	}
}

ILOCUTCALLBACK6(CtCallBack, NumVarMatrix, x, NumVarMatrix, xx, IloNumVarArray, z, NumMatrix, D, IloNum, tau, NumMatrix, pi)
{
	try
	{
	    IloEnv env = getEnv();

	    IloInt i,j,k,l,t,a,b, k1, k2, T;

	    /////////// Only add cuts before branching ///////////////////
	   if(getNnodes()==0 && nbIteration<50-eps)
	   {
		nbIteration++;

		//////////////////////////////////////////////////////////////////////////////////
		//////////// add mixing cuts at non-anticipative node ////////////////////////////
		//////////////////////////////////////////////////////////////////////////////////
		
		NumMatrix LB(env, n+1);

		for(t=1; t<=n; t++)
		{
			LB[t]=IloNumArray(env, pow(2,t-1)+1);			// Lower bound of xx[t][i]

			for(j=1; j<=pow(2,t-1); j++)
			{
				if(D[t][2*j-1]>=D[t][2*j])
				{
					if(pi[t][2*j-1]>1-tau)
					{
						if (getValue(xx[t][j])<D[t][2*j-1]-eps)
						{
							add(xx[t][j]>=D[t][2*j-1]);
							nbMix++;
						}
				
						LB[t][j]=D[t][2*j-1];

					}

					else if (pi[t][2*j-1]<=1-tau && pi[t][2*j-1]+pi[t][2*j]>1-tau)
					{
						if (getValue(xx[t][j])<D[t][2*j]-eps)
						{
							add(xx[t][j]>=D[t][2*j]);
							nbMix++;
						}

				
						IloNum ll=m/pow(2,t);

						for(k=1; k<=ll; k++)
						{
							if (getValue(xx[t][j])+(D[t][2*j-1]-D[t][2*j])*getValue(z[(2*j-2)*ll+k])<D[t][2*j-1]-eps)
							{
								add(xx[t][j]+(D[t][2*j-1]-D[t][2*j])*z[(2*j-2)*ll+k]>=D[t][2*j-1]);
								nbMix++;
							}
						}

						LB[t][j]=D[t][2*j];
					}

					else if (pi[t][2*j-1]+pi[t][2*j]<=1-tau)
					{
						IloNum ll=m/pow(2,t);

						for(k=1; k<=ll; k++)
						{
							for(i=1; i<=ll; i++)
							{
								if (getValue(xx[t][j])+(D[t][2*j-1]-D[t][2*j])*getValue(z[(2*j-2)*ll+k])+D[t][2*j]*getValue(z[(2*j-1)*ll+i])<D[t][2*j-1]-0.001)
								{
									add(xx[t][j]+(D[t][2*j-1]-D[t][2*j])*z[(2*j-2)*ll+k]+D[t][2*j]*z[(2*j-1)*ll+i]>=D[t][2*j-1]);
									nbMix++;
								}
							}
						}

						LB[t][j]=0;
					}
				}

				else if(D[t][2*j-1]<D[t][2*j])
				{
					if(pi[t][2*j]>=1-tau)
					{
						if (getValue(xx[t][j])<D[t][2*j]-eps)
						{
							add(xx[t][j]>=D[t][2*j]);
							nbMix++;
						}

						LB[t][j]=D[t][2*j];

					}

					else if (pi[t][2*j]<=1-tau && pi[t][2*j]+pi[t][2*j-1]>=1-tau)
					{
						if (getValue(xx[t][j])<D[t][2*j-1]-eps)
						{
							add(xx[t][j]>=D[t][2*j-1]);
							nbMix++;
						}

						IloNum ll=m/pow(2,t);

						for(k=1; k<=ll; k++)
						{
							if (getValue(xx[t][j])+(D[t][2*j]-D[t][2*j-1])*getValue(z[(2*j-1)*ll+k])<D[t][2*j]-eps)
							{
								add(xx[t][j]+(D[t][2*j]-D[t][2*j-1])*z[(2*j-1)*ll+k]>=D[t][2*j]);
								nbMix++;
							}
						}

						LB[t][j]=D[t][2*j-1];
					}

					else if (pi[t][2*j]+pi[t][2*j-1]<=1-tau)
					{
						IloNum ll=m/pow(2,t);

						for(k=1; k<=ll; k++)
						{
							for(i=1; i<=ll; i++)
							{
								
								if (getValue(xx[t][j])+(D[t][2*j]-D[t][2*j-1])*getValue(z[(2*j-1)*ll+k])+D[t][2*j-1]*getValue(z[(2*j-2)*ll+i])<D[t][2*j]-eps)
								{
									add(xx[t][j]+(D[t][2*j]-D[t][2*j-1])*z[(2*j-1)*ll+k]+D[t][2*j-1]*z[(2*j-2)*ll+i]>=D[t][2*j]);
									nbMix++;
								}

							}
						}

						LB[t][j]=0;

					}
				}

			}
		}


	     } 
	
	}	
	
	catch (IloException &e)
	{
		cerr<<"Error in callback:"<<e<<endl;
	}
}



int main(int argc, char **argv) 
{
   IloEnv env;
   try
   { 
	IloInt i,j,k,l,t;

	n=atoi(argv[1]);			// number of periods
	m=pow(2,n);				// number of scenarios
	IloNum tau=0.01*atoi(argv[2]);  	// service level	
	IloInt ratio = atoi(argv[3]);		// fixed to variable production cost ratio
	IloInt seed = atoi(argv[4]);		// random seed for demand and cost

	srand (seed);	

	env.out()<<"n="<<n<<" m="<<m<<" tau="<<tau<<" seed="<<seed<<" f/c="<<ratio<<endl;
	

	///////////////////// Setting up parameters and generating data ////////////////////////////

	NumMatrix d(env, n+1);		// individual period demand
	NumMatrix c(env, n+1);         	// variable production cost
	NumMatrix f(env, n+1);    	// fixed production cost
	NumMatrix h(env, n+1);		// inventory holding cost
	NumMatrix D(env, n+1);	        // cumulative demand
	NumMatrix pi(env, n+1);		// probability for each scenario


	for (j=0; j<=n; j++)
	{
		d[j] = IloNumArray (env, pow(2,j)+1);
		c[j] = IloNumArray (env, pow(2,j)+1);
		f[j] = IloNumArray (env, pow(2,j)+1);
		h[j] = IloNumArray (env, pow(2,j)+1);
		D[j] = IloNumArray (env, pow(2,j)+1);
		pi[j] = IloNumArray (env, pow(2,j)+1);
	}

	pi[0][1]=1;
	D[0][1]=0;

	for (j=1; j<=n; j++)
	{
		for (i=1; i<=pow(2,j-1); i++)
		{
			IloNum P=rand () %200 + 1;
			P=P/200;
			pi[j][2*i-1]=pi[j-1][i]*P;
			pi[j][2*i]=pi[j-1][i]*(1-P);
		}
	}
	
	for (j=1; j<=n; j++)
	{
		for (i=1; i<=pow(2,j-1); i++)
		{
			d[j][2*i-1]=rand() %50+50;
			d[j][2*i]=rand() % 50+50;

			c[j][2*i-1]=rand() %10+10;	
			c[j][2*i]=rand() % 10+10;

			h[j][2*i-1]=0.1*c[j][2*i-1];
			h[j][2*i]=0.1*c[j][2*i];

			f[j][2*i-1]=ratio*c[j][2*i-1];
			f[j][2*i]=ratio*c[j][2*i];
			
			D[j][2*i-1]=D[j-1][i]+d[j][2*i-1];
			D[j][2*i]=D[j-1][i]+d[j][2*i];

		}
	}

	
	IloNum MaxD=D[n][m];

	for(i=1; i<=m-1; i++)
	{
		if(MaxD<D[n][i]-eps)
		{
			MaxD=D[n][i];
		}
	}


	for (t=1; t<=n; t++)
	{	
		for (j=1; j<=pow(2,t); j++)
		{
			env.out()<<D[t][j]<<" ";
		}
		env.out()<<endl;
	}
	env.out()<<endl;

	/////////////// Setting up variables ////////////////////////////
	
	NumVarMatrix x(env, n+1);			// production variable
	NumVarMatrix xx(env, n+1);	
        NumVarMatrix s(env, n+1);			// inventory variable
	NumVarMatrix y(env, n+1);			// set up variable, ==1 if set up at time j in scenario k
    	IloNumVarArray z(env, m+1, 0, 1, ILOINT);	//stock out variable, ==1 if stockout in scenario k

	for (j=1; j<=n; j++)
	{
		x[j]= IloNumVarArray (env, pow(2,j-1)+1, 0, IloInfinity, ILOFLOAT);
		xx[j]= IloNumVarArray (env, pow(2,j-1)+1, 0, IloInfinity, ILOFLOAT);
	  	s[j]= IloNumVarArray (env, pow(2,j)+1, 0, IloInfinity, ILOFLOAT);
	  	y[j]= IloNumVarArray (env, pow(2,j-1)+1, 0, 1, ILOINT);
	}


	/////////////////// Build model ///////////////////////////////////

	IloModel model(env); 
	IloCplex cplex(env);

	///////////////////// Constraints //////////////////////////////////
	

	model.add(xx[1][1]==x[1][1]);

	for(j=2; j<=n; j++)
	{
		for(i=1; i<=pow(2,j-2); i++)
		{
			model.add(xx[j][2*i-1]==xx[j-1][i]+x[j][2*i-1]);
			model.add(xx[j][2*i]==xx[j-1][i]+x[j][2*i]);
		}
	}
	
	// knapsack constraint
	IloExpr lhs(env);
	for (k=1; k<=m; k++)
	{
		lhs += pi[n][k]*z[k];
	}
	model.add (lhs <= 1-tau);
	lhs.end();

	// Probability constraint
	for (j=1; j<=n; j++)
	{
	 	for (i=1; i<=pow(2,j-1); i++)
	 	{
			IloInt length=m/pow(2,j);
			for(k=1; k<=length; k++)
			{
				model.add (xx[j][i]+D[j][2*i-1]*z[(2*i-2)*length+k] >= D[j][2*i-1]);
				model.add (xx[j][i]+D[j][2*i]*z[(i*2-1)*length+k] >= D[j][2*i]);


			}
		}
	}


	// endhorizon constraint
	
	for(k=1; k<=pow(2,n-1); k++)
	{
		model.add(xx[n][k]>=D[n][2*k-1]);
		model.add(xx[n][k]>=D[n][2*k]);
	}

	// Inventory constraint {j in 1..t, k in 1..s}: i[j,k] >= sum {l in 1..j} (x[l,k]-demand[l,k]);

	for (j=1; j<=n; j++)
	{
		for (k=1; k<=pow(2,j-1); k++)
	  	{
			model.add (s[j][2*k-1] >= xx[j][k] - D[j][2*k-1]);
			model.add (s[j][2*k] >= xx[j][k] - D[j][2*k]);
		}
		
    	}
	// fixed cost constraint
	for(j=1; j<=n; j++)
	{
		for(k=1; k<=pow(2,j-1); k++)
		{
			model.add(x[j][k] <= (MaxD-20)*y[j][k]);
		}
	}

	/////////////// Objective Function: Minimize Cost //////////////////

	IloExpr obj(env);	

 	for (j=1; j<=n; j++)
	{
		for(k=1; k<=pow(2, j-1); k++)
		{
			obj += pi[j][2*k-1]*(c[j][2*k-1]*x[j][k]+h[j][2*k-1]*s[j][2*k-1]+f[j][2*k-1]*y[j][k]);
			obj += pi[j][2*k]*(c[j][2*k]*x[j][k]+h[j][2*k]*s[j][2*k]+f[j][2*k]*y[j][k]);
		}
	}
        model.add(IloMinimize(env, obj));
	obj.end();

	////////////////// Build model & solve ////////////////////////////////

        env.out() << "constructing cut callback ..." << endl;
      	env.out() << "extracting model ..." << endl;
	cplex.extract(model);

	cplex.setParam(IloCplex::TiLim, 3600);  
	cplex.setParam(IloCplex::EpGap, 0);
	cplex.setParam(IloCplex::MIPSearch, 1);  
	cplex.setParam(IloCplex::Threads, 1);  
//	cplex.setParam(IloCplex::HeurFreq, -1);
	cplex.setParam(IloCplex::PreInd,0);
	cplex.setParam(IloCplex::MIPDisplay, 4);  

	cplex.use(CtCallBack(env, x, xx, z, D, tau, pi));


	

	///////////////// Close non-structural cuts (non-flow) ///////////////////////////

//	cplex.setParam(IloCplex::FracCuts,-1);		cplex.setParam(IloCplex::DisjCuts,-1); 		cplex.setParam(IloCplex::ZeroHalfCuts,-1);
//	cplex.setParam(IloCplex::MIRCuts,-1); 		cplex.setParam(IloCplex::ImplBd,-1);		cplex.setParam(IloCplex::GUBCovers,-1); 
//	cplex.setParam(IloCplex::Cliques,-1);		cplex.setParam(IloCplex::Covers,-1); 
	
 	
	///////////////// Open non-structural cuts (non-flow) ///////////////////////////

	cplex.setParam(IloCplex::FracCuts,2);		cplex.setParam(IloCplex::DisjCuts,2); 		cplex.setParam(IloCplex::ZeroHalfCuts,2);
	cplex.setParam(IloCplex::MIRCuts,2); 		cplex.setParam(IloCplex::ImplBd,2);		cplex.setParam(IloCplex::GUBCovers,2); 
	cplex.setParam(IloCplex::Cliques,2);		cplex.setParam(IloCplex::Covers,2); 
		
	///////////////// Close structural cuts (flow related) ///////////////////////////
	    	     	
//	cplex.setParam(IloCplex::FlowCovers,-1);     	cplex.setParam(IloCplex::FlowPaths,-1);	 		
//	cplex.setParam(IloCplex::MCFCuts, -1);

	///////////////// Open structural cuts (flow related) ///////////////////////////
		       	
	cplex.setParam(IloCplex::FlowCovers,2);     	cplex.setParam(IloCplex::FlowPaths,2);
	cplex.setParam(IloCplex::MCFCuts, 2);    		   	
    			
    			



        clock_t start, end;
     	double cpu_time_used;

	start = clock();
	cplex.solve();
	env.out() << "solving model ...\n";
	end = clock();

	cpu_time_used = ((double) (end - start)) / CLOCKS_PER_SEC;

	if (cplex.getStatus() == IloAlgorithm::Infeasible)
                         env.out() << "No Solution" << endl;

	
	//////////////////// Print results /////////////////////////////////////
	env.out()<<endl<<"Mixing cuts applied: "<<nbMix<<endl;
	env.out()<<"Continuous mixing cuts applied: "<<nbCmix<<endl;

	IloNum ser=0;

	for(i=1; i<=m; i++)
	{	
		ser=ser+(1-cplex.getValue(z[i]))*pi[n][i];
	}
	
	
	env.out() << endl << "Minimum cost is: " << cplex.getObjValue() << endl;

	env.out() << "Solution status = " << cplex.getStatus() << endl;
     		
	env.out() << "End gap  = " << 100*(cplex.getObjValue()-cplex.getBestObjValue())/cplex.getObjValue()<<"%"<< endl;
	
	env.out() <<"# of nodes: "<<cplex.getNnodes() << endl;
	
	env.out()<<"time="<<cpu_time_used<<endl;


	env.out()<<endl;

	for(t=1; t<=n; t++)
	{	
		for(j=1; j<=pow(2,t-1); j++)
		{
       			env.out()<<cplex.getValue(x[t][j])<<", ";
		}
		env.out()<<endl;

		for(j=1; j<=pow(2,t-1); j++)
		{
       			env.out()<<cplex.getValue(xx[t][j])<<", ";
		}
		env.out()<<endl;
	
		for (j=1; j<=pow(2,t-1); j++)
              	{
			env.out()<<cplex.getValue(y[t][j])<<", ";
		}
		env.out()<<endl;
				
		for (j=1; j<=pow(2,t); j++)
              	{
			env.out()<<cplex.getValue(s[t][j])<<", ";
		}
		env.out()<<endl<<endl;
	}

    }

    catch (IloException& ex) 
    {
    	cerr << "Error: " << ex << endl;
    }
    catch (...) 
    {
   	cerr << "Error" << endl;
    }

    env.end();
    return 0;
}

