#define EPS 0.0000001

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

ILOSTLBEGIN
typedef IloArray<IloNumArray> NumArray2;

IloInt numCuts = 0;

ILOCUTCALLBACK3(lotSizingCuts, NumArray2, D, IloNumVarArray, product, IloNumVarArray, setup) {
     try{
        IloEnv env = getEnv();
        IloInt NT = product.getSize();

        for (int t = 0; t < NT; t++){
            IloExpr lhs(env);
            IloNum value = 0;
            for (int s = 0; s <= t; s++){
                if (getValue(product[s]) > D[s][t]*getValue(setup[s]) + EPS){
                    lhs += D[s][t] * setup[s];
                    value += D[s][t] * getValue(setup[s]);
                }
                else{
                    lhs += product[s];
                    value += getValue(product[s]);
                }
            }
            if (value < D[0][t] - EPS){
                add(lhs >= D[0][t]);
                numCuts ++;

	        env.out()<<lhs<< " >= "<<D[0][t]<<endl; 

            }
            lhs.end();
        }
    }

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

}
int main(int argc, char **argv) {
   IloEnv env;
   try {
     IloInt s,t,k;

     IloInt NT = 6;
     IloNumArray DEMAND(env, NT, 1, 3, 5, 3, 4, 2);
     IloNumArray SETUPCOST(env, NT, 17, 16, 11, 6, 9, 6);
     IloNumArray PRODCOST(env, NT, 5, 3, 2, 1, 3, 1);

     IloNumVarArray product(env, NT, 0, IloInfinity, ILOFLOAT);
     IloNumVarArray setup(env, NT, 0, 1, ILOINT);
     NumArray2 D(env, NT);
     IloModel model(env);

     for(s = 0; s < NT; s++){
       D[s] = IloNumArray(env, NT);
       for(t = 0; t < NT; t++)
         for(k = s; k <= t; k++)
           D[s][t] += DEMAND[k];
     }

     IloExpr obj(env);
     obj = IloScalProd(SETUPCOST,setup) + IloScalProd(PRODCOST,product);
     model.add(IloMinimize(env, obj));

     for(t = 0; t < NT; t++){
       IloExpr lhs(env), rhs(env);
       for( s = 0; s <= t; s++){
         lhs += product[s];
         rhs += DEMAND[s];
       }
       model.add( lhs >= rhs );
       lhs.end(); rhs.end();
     }
     for(t = 0; t < NT; t++)
       model.add( product[t] <= D[t][NT-1]*setup[t] );


     IloCplex cplex(model);

     cplex.use(lotSizingCuts(env,D,product,setup));

     cplex.exportModel("ls-model.lp");

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

     ofstream logfile("cplex.log");
     cplex.setOut(logfile);
     cplex.setWarning(logfile);


        
     cplex.setParam(IloCplex::TiLim, 1800);
     cplex.setParam(IloCplex::TreLim, 100);
     cplex.setParam(IloCplex::NodeLim, 1000000);

     cplex.setParam(IloCplex::MIPDisplay, 2);
     cplex.setParam(IloCplex::HeurFreq, -1);
     cplex.setParam(IloCplex::PreInd,0);
     cplex.setParam(IloCplex::Covers,-1);
     cplex.setParam(IloCplex::Cliques,-1);
     cplex.setParam(IloCplex::FlowCovers,-1);
     cplex.setParam(IloCplex::GUBCovers,-1);
     cplex.setParam(IloCplex::FracCuts,-1);
     cplex.setParam(IloCplex::MIRCuts,-1);
     cplex.setParam(IloCplex::FlowPaths,-1);
     cplex.setParam(IloCplex::ImplBd,-1);
     cplex.setParam(IloCplex::DisjCuts,-1);




     cplex.solve();


     env.out() << endl << "Optimal value: " << cplex.getObjValue() << endl;
     env.out()<< "Number of Cuts added: " << numCuts <<endl;
     env.out()<< "Number of B&B nodes: " <<cplex.getNnodes() << endl;  
     for(t = 0; t < NT; t++)
       if (cplex.getValue(setup[t]) >= 1 - eps)
            env.out() << "At time " << t << " produced " <<
            cplex.getValue(product[t]) << " " << 
            cplex.getValue(setup[t]) << endl;
   }
   catch(IloException& e) {
     cerr  << " ERROR: " << e << endl;
   }
   catch(...) {
     cerr  << " ERROR" << endl;
   }
   env.end();
   return 0;
}