// This is the last version of the code mainmpi.cpp. All rights reserved.
// Copyright (C) 2012 by U. Aldasoro, M.A. Garin, M. Merino and G. Perez,
// from UPV/EHU.
// No part of this code may be reproduced, modified or transmitted, in any
// form or by any means without the prior written permission of the authors. 

// Cluster models are in mps format, as external files

#include "itoa.h"
#include <mpi.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <stdio.h>
#include <iostream>
#include <fstream>
using namespace std;
#include "ClpSimplex.hpp"
#include "CoinHelperFunctions.hpp"
#include "CoinBuild.hpp"
#include "CoinModel.hpp"
#include <iomanip>
#include <cassert>
#include "OsiClpSolverInterface.hpp"
#include "CbcModel.hpp"

#define rmaxmin (1)		// (1): Maximum; (-1): Minimum
#define nmodel 44		// Number of cluster models
#define useCplex 1		// Select solver.	0: COIN-OR; 1: Cplex
#define threadsCplex 2  // Number of parallel threads used by each MPI thread

#ifdef useCplex
	#include "OsiCpxSolverInterface.hpp"
#endif

int main(int argc, char **argv)
{    
	// STEP 0 - DECLARING OPTIMIZATION AND MPI VARIABLES
	int imod,i,j,loc,pid,npr,lp,original_rank,original_size,error,ncols;
	int assignment[nmodel],inicial[nmodel],nonzero[nmodel],numres[nmodel],
		numvar[nmodel],numvarint[nmodel],assignmentX0[nmodel],iniX0[nmodel],
		istruef969[nmodel];
	double zq[nmodel];

	MPI_Group  orig_group,new_group;		//MPI groups
	MPI_Comm   new_comm;					//MPI communicator


	// STEP 1 - DEFINITION OF THE GLOBAL ENVIRONMENT
	MPI_Init(&argc,&argv); // Beginning of the MPI environment
		MPI_Comm_size(MPI_COMM_WORLD, &original_size); // Total number of threads
		MPI_Comm_rank(MPI_COMM_WORLD, &original_rank); // Who am I?

		// 1.1 - Creating a new communicator considering active threads

		// 1.1.1- Basic variables
		int ranks1[nmodel];
		int *ranks2;  ranks2=new int[original_size];
		for(i=0;i<nmodel;i++) {ranks1[i]=i;}
		for(i=0;i<original_size;i++) {ranks2[i]=i;}

		// 1.1.2- Extract handle of global group from MPI_COMM_WORLD using MPI_Comm_group
		MPI_Comm_group(MPI_COMM_WORLD, &orig_group);

		// 1.1.3- Form new group as a subset of global group using MPI_Group_incl
		if (nmodel < original_size){MPI_Group_incl(orig_group, nmodel, ranks1, &new_group); 
		} else {MPI_Group_incl(orig_group, original_size, ranks2, &new_group);}

		// 1.1.4- Create new communicator for new group using MPI_Comm_create
		MPI_Comm_create(MPI_COMM_WORLD, new_group, &new_comm);

		// 1.1.5- Determine new rank in new communicator using MPI_Comm_rank
		MPI_Group_rank (new_group, &pid);						
		MPI_Group_size (new_group, &npr);

		// 1.1.6 - Creating the output file
		ofstream solution("CPLEX_COIN_MPI_solution.dat",ios::out); //Output file

		// 1.2 - Only active threads follow the solving process

		if (original_rank < npr) {

			// 1.2.1 - Assigning the work load
			for(i=0;i<npr;i++) assignment[i]=int(nmodel/npr);
			for(i=0;i<(nmodel-npr*int(nmodel/npr));i++) assignment[i]=assignment[i]+1;

			inicial[0]=0;
			for(i=1;i<npr;i++) inicial[i]=inicial[i-1]+assignment[i-1];

			// 1.2.2 - Declaring local variables
			int ncols_max_loc=0;
			double *zq_loc;  zq_loc=new double[assignment[pid]];
			int *nonzero_loc;  nonzero_loc=new int[assignment[pid]];
			int *numres_loc;  numres_loc=new int[assignment[pid]];
			int *numvar_loc;  numvar_loc=new int[assignment[pid]];
			int *numvarint_loc;  numvarint_loc=new int[assignment[pid]];
			int *istruef969_loc;  istruef969_loc=new int[assignment[pid]];	

			// 1.2.3 - Creating the output file
			if (pid == 0) {
				if (useCplex == 0) {solution<<"*** SOLVER: COIN-OR ***"<<"\n\n";} 
				else {solution<<"*** SOLVER: CPLEX *** \n\n  +++ NUMBER OF " << 
					"PARALLEL THREADS ON CPLEX = "<<threadsCplex<<" +++ \n\n";}

				solution<<"Beginning of CPLEX_COIN_MPI_solution.dat (output of mainmpi.cpp) \n  Number of threads = "<<original_size<<"\n Number of submodels = "<<nmodel<<"\n Number of active threads = "<<npr<<"\n";
				for(i=0;i<npr;i++) {
					solution<<" Number of submodels solved by thread "<<i<<" = "<< 
						assignment[i]<<", starting by submodel = "<< inicial[i]+1<<"\n";}
			}

			// 1.2.4 - Creating the COIN-OR or CPLEX/COIN-OR models
			OsiClpSolverInterface *sol0;
			sol0=new OsiClpSolverInterface[assignment[pid]];

			OsiCpxSolverInterface *sol1;
			sol1=new OsiCpxSolverInterface[assignment[pid]];
			for(loc=0;loc<assignment[pid];loc++){
				CPXENVptr env = sol1[loc].getEnvironmentPtr();
				CPXsetintparam(env, CPX_PARAM_THREADS, threadsCplex);
			}

			for(loc=0;loc<assignment[pid];loc++) {

				imod=inicial[pid]+loc;			const char* final;
				final="";						char buffer[33];
				final=itoa(imod+1, buffer,10); //convert in char
				char model[80];					strcpy (model,"Cluster");
				strcat (model,final);			puts (model);

				if (useCplex == 0) {
					sol0[loc].setObjSense(rmaxmin);
					sol0[loc].readMps(model); //Read cluster submodel
					nonzero_loc[loc] = sol0[loc].getNumElements();
					numres_loc[loc] = sol0[loc].getNumRows();
					numvar_loc[loc] = sol0[loc].getNumCols(); 
					for(j=0;j<sol0[loc].getNumCols();j++) 
						if(sol0[loc].isInteger(j)) numvarint_loc[loc]++;
					if(sol0[loc].getNumCols() > ncols_max_loc)
						{ncols_max_loc = sol0[loc].getNumCols();}
				} else {
					sol1[loc].setObjSense(rmaxmin);
					sol1[loc].readMps(model); //Read cluster submodel
					nonzero_loc[loc] = sol1[loc].getNumElements();
					numres_loc[loc] = sol1[loc].getNumRows();
					numvar_loc[loc] = sol1[loc].getNumCols(); 
					for(j=0;j<sol1[loc].getNumCols();j++) 
						if(sol1[loc].isInteger(j)) numvarint_loc[loc]++;
					if(sol1[loc].getNumCols() > ncols_max_loc)
						{ncols_max_loc = sol1[loc].getNumCols();}
				}
			}

			if (npr > 1) {
				MPI_Allreduce(&ncols_max_loc,&ncols, 1, MPI_INT, MPI_MAX,new_comm);
			} else {
				ncols=ncols_max_loc;}

			double **x0;x0=new double*[nmodel];     
            for(i=0;i<nmodel;i++) x0[i]=new double [ncols]; 

			double *x0_loc;  x0_loc=new double[assignment[pid]*ncols];
			double *x0vector;  x0vector=new double[nmodel*ncols];

			for(imod=0;imod<assignment[pid];imod++) {
				for(j=0;j<nmodel;j++) {
						x0_loc[imod*ncols+j]=0;
				}
			}

			for(i=0;i<npr;i++) {
				iniX0[i]=inicial[i]*ncols; assignmentX0[i]=assignment[i]*ncols;}

			// STEP 2 - PRESOLVE ASSIGNED MODELS
			error=0;

			for(loc=0;loc<assignment[pid];loc++)
			{
				if (useCplex == 0) {
					imod=inicial[pid]+loc;			
					CbcModel pm0(sol0[loc]);
					istruef969_loc[loc]=0;			sol0[loc].initialSolve();

					if(sol0[loc].isProvenPrimalInfeasible()) {zq_loc[loc]=rmaxmin*1.e25; 
					istruef969_loc[loc]=1;}

					if(istruef969_loc[loc]==0){
						if(sol0[loc].isProvenDualInfeasible() ){
							zq_loc[loc]=-1.0*rmaxmin*1.e25;	istruef969_loc[loc]=2;}

						if(!sol0[loc].isProvenOptimal()) {
							zq_loc[loc]=-1.0*rmaxmin*1.e25;	istruef969_loc[loc]=2;}
					}
				} else {
					imod=inicial[pid]+loc;			
					istruef969_loc[loc]=0;			sol1[loc].initialSolve();

					if(sol1[loc].isProvenPrimalInfeasible()) {zq_loc[loc]=rmaxmin*1.e25; 
					istruef969_loc[loc]=1;}

					if(istruef969_loc[loc]==0){
						if(sol1[loc].isProvenDualInfeasible() ){
							zq_loc[loc]=-1.0*rmaxmin*1.e25;	istruef969_loc[loc]=2;}

						if(!sol1[loc].isProvenOptimal()) {
							zq_loc[loc]=-1.0*rmaxmin*1.e25;	istruef969_loc[loc]=2;}
					}
				}
			}

			// STEP 3 - GLOBAL MPI COMMUNICATION (GATHERING PRESOLVE INFORMATION)
			if (npr > 1) {
				MPI_Allgatherv(&istruef969_loc[0], assignment[pid], MPI_INT,
					&istruef969[0],assignment, inicial, MPI_INT, new_comm);
				MPI_Gatherv(&zq_loc[0], assignment[pid], MPI_INT, &zq[0],
					assignment, inicial, MPI_INT, 0, new_comm);
			} else {
				for(imod=0;imod<nmodel;imod++) {
					istruef969[imod]=istruef969_loc[imod];		zq[imod]=zq_loc[imod];
				}
			}	

			for(imod=0;imod<nmodel;imod++) {
				if(istruef969[imod] == 1){
					solution<<"\n Cluster submodel "<<imod+1 
						<<" is primal infeasible. ObjValue = "<<zq[imod];	error=1;}
				if(istruef969[imod] == 2){
					solution<<"\n Cluster submodel "<<imod+1
						<<" is unbounded. ObjValue = "<<zq[imod];	error=1;}
			}

			// STEP 4 - CHECK PRESOLVE STATUS
			if(error == 0){

				// STEP 5- SOLVE ASSIGNED MODELS
				if (useCplex == 0) {
					for(loc=0;loc<assignment[pid];loc++)
					{
						imod=inicial[pid]+loc;		CbcModel pm0(sol1[loc]);
						pm0.branchAndBound();

						if(!pm0.isProvenOptimal()) {
							zq_loc[loc]=-1.0*rmaxmin*1.e25;		istruef969_loc[loc]=3;}

						if(fabs(pm0.getBestPossibleObjValue())>1.e24){
							zq_loc[loc]=-1.0*rmaxmin*1.e25;		istruef969_loc[loc]=4;}

						if(istruef969_loc[loc] == 0) {
							zq_loc[loc]=pm0.getObjValue();
							for(j=0;j<pm0.getNumCols();j++) {
								x0_loc[loc*ncols+j]=pm0.getColSolution()[j];}}
					}
				} else {
					for(loc=0;loc<assignment[pid];loc++)
					{
						imod=inicial[pid]+loc;		sol1[loc].branchAndBound();

						if(!sol1[loc].isProvenOptimal()) {
							zq_loc[loc]=-1.0*rmaxmin*1.e25;	istruef969_loc[loc]=3;}

						if(fabs(sol1[loc].getObjValue())>1.e24){
							zq_loc[loc]=-1.0*rmaxmin*1.e25;	istruef969_loc[loc]=4;}

						if(istruef969_loc[loc] == 0) {
							zq_loc[loc]=sol1[loc].getObjValue();
							for(j=0;j<sol1[loc].getNumCols();j++) 
								x0_loc[loc*ncols+j]=sol1[loc].getColSolution()[j];}
					}
				}

				// STEP 6- GLOBAL MPI COMMUNICATION (GATHERING SOLVE INFORMATION)
				if (npr > 1) {
					MPI_Gatherv(&zq_loc[0],assignment[pid],MPI_DOUBLE,&zq[0],
						assignment,inicial,MPI_DOUBLE,0,new_comm);
					MPI_Gatherv(&nonzero_loc[0],assignment[pid],MPI_INT,&nonzero[0],
						assignment,inicial,MPI_INT,0,new_comm);
					MPI_Gatherv(&numvar_loc[0],assignment[pid],MPI_INT,&numvar[0],
						assignment,inicial,MPI_INT,0,new_comm);
					MPI_Gatherv(&numvarint_loc[0],assignment[pid],MPI_INT,&numvarint[0],
						assignment,inicial,MPI_INT,0,new_comm);
					MPI_Gatherv(&numres_loc[0],assignment[pid],MPI_INT,&numres[0],
						assignment,inicial,MPI_INT,0,new_comm);
					MPI_Gatherv(x0_loc,assignmentX0[pid],MPI_DOUBLE,x0vector,
						assignmentX0,iniX0,MPI_DOUBLE,0,new_comm);
					MPI_Gatherv(&istruef969_loc[0],assignment[pid],MPI_INT,&istruef969[0],
						assignment,inicial,MPI_INT,0,new_comm);
				} else {
					for(imod=0;imod<nmodel;imod++) {
						zq[imod]=zq_loc[imod];				nonzero[imod]=nonzero_loc[imod];		
						numvar[imod]=numvar_loc[imod];		numres[imod]=numres_loc[imod];
						istruef969[imod]=istruef969_loc[imod];
						for(j=0;j<numvar[imod];j++){
							x0vector[imod*ncols+j]=x0_loc[imod*ncols+j];
						}
					}
				}

				for(imod=0;imod<nmodel;imod++) {
					if((istruef969[imod] == 3) | (istruef969[imod] == 4)){
						solution<<"\n Cluster submodel "<<imod+1 <<" is unbounded ";}
				}
			}

			// STEP 7 - PRINT RESULTS (if pid == 0)
			if (pid == 0) {

				solution<<"\n\n [RESULTS 1:] PROPERTIES OF SOLVED MODELS ";
				for(imod=0;imod<nmodel;imod++){
					solution<<"\n \n Cluster submodel "<<imod+1<<" of "<<nmodel;        
					solution<< "\n Number of variables: "; solution<< numvar[imod];   
					solution<< "\n Number of constraints: "; solution<< numres[imod];
					solution<< "\n Number of nonzero elements: "<<nonzero[imod];
					for(j=0;j<numvar[imod];j++){
						x0[imod][j]=x0vector[imod*ncols+j];
						//solution<< "\n x[ "<<imod<<","<<j<<"]"<<x0[imod][j]<<" \n"; // To print x0
					}
				}

				solution<<"\n\n\n [RESULTS 2:] OBJECTIVE FUNCTION VALUES \n";

				for(imod=0;imod<nmodel;imod++) {
					solution<<"\n Optimal objective value Cluster submodel "<<imod+1
						<<" of "<<nmodel <<" is "<<zq[imod];}   
			}

		}

	// STEP 8 - EXECUTION ENDS IN ALL PROCESSORS
	solution.close();
	MPI_Finalize(); //End of the MPI environment
	return 0;
}