MERA_CPLEX.txt

/*
 * This program is to solve the Most Economic Redudent Assignment problem.
 * Author: Dr. Haibin Zhu, 2018
 * Please cite:
 [1] H. Zhu, “The Most Economical Redundant Assignment”, The IEEE Int'l Conf. on Systems, Man and Cybernetics (SMC'19), Bari, Italy, Oct. 6-9, 2019, pp. 146-151.
 [2] H. Zhu, E-CARGO and Role-Based Collaboration: Modeling and Solving Problems in the Complex World, Wiley-IEEE Press, NJ, USA, Dec. 2021. 
 [3] H. Zhu, M.C. Zhou, and R. Alkins, “Group Role Assignment via a Kuhn-Munkres Algorithm-based Solution”, IEEE Trans. on Systems, Man, and Cybernetics, Part A: Systems and Humans, vol. 42, no. 3, May 2012, pp. 739-750.
 [4]	H. Zhu, and M. Zhou, “Role-Based Collaboration and its Kernel Mechanisms,” IEEE Trans. on Systems, Man, and Cybernetics, Part C: Applications and Reviews, vol. 36, no. 4, July. 2006, pp. 578-589.
 * 
*/

import ilog.concert.*;
import ilog.cplex.*;

import java.io.BufferedWriter;
import java.io.IOException;
import java.text.DecimalFormat;
import java.util.*;

class GMEO_ILOG {
	
	private int m;	//number of agents
	private int n;	//number of roles
	private int limit;	//number of training roles
	private int[] L;	//Requirement array
	private int[][] A;	//Requirement array
	private double [] QV;	//Requirement array
	
	
	DecimalFormat df = new DecimalFormat("0.00");
	
	double optimized_result = 0;
	boolean bILOG_result;
	
	public GMEO_ILOG(int nagent, int nrole, double [][] QD, int[]RA)
	{
		m = nagent;
		n = nrole;
				
		L = new int[n];
		L = RA;
		System.out.println("Requirement Array: ");
		for(int i=0; i<n; i++)		System.out.print(L[i]+"	");
		System.out.print("\n");	
		A =new int [m][n];
		QV = new double [m*n];
		for(int i=0; i<m; i++)	for(int j=0; j<n; j++) QV[n*i+j]=QD[i][j];		
	}
	
	public boolean resolve(int[][]TR, int lmt)
	{
		try
		{
			limit = lmt;
			//Creat cplex obj
			IloCplex cplex = new IloCplex();	//initialize the cplex object
			
			IloIntVar[]x = cplex.intVarArray(m*n, 0, 1);	//initialize the variables array under cplex.
			IloIntVar[]y = cplex.intVarArray(m*(m-1)*n, 0, 1);	//initialize the variables array under cplex.
			
			cplex.addMaximize(cplex.scalProd(x,QV));	//add the optimize objective to cplex.
			
			//Add Constrains:								
			//Constrain type 2, T<=QD
			for (int a = 0; a<m; a++)
				for (int j = 0; j<n; j++)
				{	
					for (int i = 0; i<a; i++){
						IloLinearNumExpr exprReqConstrain = cplex.linearNumExpr();
						exprReqConstrain.addTerm(1, y[a*(m-1)*n+j+i*n]);
						exprReqConstrain.addTerm(-1, x[j+i*n]);
						cplex.addLe(exprReqConstrain,0);
					}
					for (int i = a+1; i<m; i++){
						IloLinearNumExpr exprReqConstrain = cplex.linearNumExpr();
						exprReqConstrain.addTerm(1, y[a*(m-1)*n+j+(i-1)*n]);
						exprReqConstrain.addTerm(-1, x[j+i*n]);
						cplex.addLe(exprReqConstrain,0);
					}
				}
			//Constrain type 3, L: Add role requirement constrains for each QD without agent a. 
			//the number of people assigned on each role should meet the requirement on that role.
			//Hence, n constrains will be added.
			for (int a = 0; a<m; a++){	
				for (int j = 0; j<n; j++)
				{	
					IloLinearNumExpr exprReqConstrain = cplex.linearNumExpr();
					for (int i = 0; i<a; i++)	exprReqConstrain.addTerm(1, y[a*(m-1)*n+j+i*n]);
					for (int i = a+1; i<m; i++)	exprReqConstrain.addTerm(1, y[a*(m-1)*n+j+(i-1)*n]);					
					cplex.addEq(exprReqConstrain, L[j]);
				}
			}
			//Constrain type 4: an agent can only be assigned with one role.
			for (int a = 0; a<m; a++)	
				for(int i=0; i<m; i++)
				{	IloLinearNumExpr exprUniConstrain = cplex.linearNumExpr();
					if (i<a) for(int j = 0; j<n; j++)	exprUniConstrain.addTerm(1, y[a*(m-1)*n+n*i+j]);
					if (i>a) for(int j = 0; j<n; j++)	exprUniConstrain.addTerm(1, y[a*(m-1)*n+n*(i-1)+j]);
					cplex.addLe(exprUniConstrain, 1);	
				}
			//Constrain type 5: The total number of assigned training roles is less than nt
			IloLinearNumExpr exprConstraint = cplex.linearNumExpr();
			for(int i=0; i<m; i++)
				for(int j = 0; j<n; j++)	exprConstraint.addTerm(1, x[n*i+j]);
			cplex.addLe(exprConstraint, limit);	

			
			System.out.println("Before Resolve! ");

			//Solve LP
			//long t1 = System.nanoTime();
			if (cplex.solve()) 
			{
				bILOG_result = true;
				optimized_result = cplex.getObjValue();
				cplex.output().println("Solution status = " + cplex.getStatus());
				cplex.output().println("Solution value = " + cplex.getObjValue());
				
				double[] val = cplex.getValues(x);
				int ncols = cplex.getNcols();
				cplex.output().println("Num COL: " + ncols);
				
				cplex.output().println("Result Table: " );
				for (int j=0; j<ncols; j++)
				{
					if (j/n==m) break;
					A[j/n][j%n] = (int)val[j]; 
					System.out.print(A[j/n][j%n] + " ");
					TR[j/n][j%n] = A[j/n][j%n];
					//System.out.print(val[j]+ "	");
					if ((j+1)%(n) == 0) {System.out.print("\n");}	
				}
//				double[] val1 = cplex.getValues(y);
//				System.out.print("Etta ="+(int)val1[0]);

				cplex.end();
				bILOG_result = true;				
			}
			else
			{
				cplex.end();
				bILOG_result = false;
			}	
			//long t2 = System.nanoTime();
			//time[0] = (t2-t1)/1000000;
		}
		catch (IloException e){System.err.println("Concert exception" + e + " caught");}
		
		return(bILOG_result);
	}
	
	public double getOptimizedResult()
	{
		return optimized_result;
		
	}
}
class RatedMunkres {

    static public int[][] computeAssignments(double[][] matrix) {


        // subtract minumum value from rows and columns to create lots of zeroes
        reduceMatrix(matrix);


        // non negative values are the index of the starred or primed zero in the row or column
        int[] starsByRow = new int[matrix.length]; Arrays.fill(starsByRow,-1);
        int[] starsByCol = new int[matrix[0].length]; Arrays.fill(starsByCol,-1);
        int[] primesByRow = new int[matrix.length]; Arrays.fill(primesByRow,-1);

        // 1s mean covered, 0s mean not covered
        int[] coveredRows = new int[matrix.length];
        int[] coveredCols = new int[matrix[0].length];

        // star any zero that has no other starred zero in the same row or column
        initStars(matrix, starsByRow, starsByCol);
        coverColumnsOfStarredZeroes(starsByCol,coveredCols);

        while (!allAreCovered(coveredCols)) {

            int[] primedZero = primeSomeUncoveredZero(matrix, primesByRow, coveredRows, coveredCols);

            while (primedZero == null) {
                // keep making more zeroes until we find something that we can prime (i.e. a zero that is uncovered)
                makeMoreZeroes(matrix,coveredRows,coveredCols);
                primedZero = primeSomeUncoveredZero(matrix, primesByRow, coveredRows, coveredCols);
            }

            // check if there is a starred zero in the primed zero's row
            int columnIndex = starsByRow[primedZero[0]];
            if (-1 == columnIndex){

                // if not, then we need to increment the zeroes and start over
                incrementSetOfStarredZeroes(primedZero, starsByRow, starsByCol, primesByRow);
                Arrays.fill(primesByRow,-1);
                Arrays.fill(coveredRows,0);
                Arrays.fill(coveredCols,0);
                coverColumnsOfStarredZeroes(starsByCol,coveredCols);
            } else {

                // cover the row of the primed zero and uncover the column of the starred zero in the same row
                coveredRows[primedZero[0]] = 1;
                coveredCols[columnIndex] = 0;
            }
        }

        // ok now we should have assigned everything
        // take the starred zeroes in each column as the correct assignments

        int[][] retval = new int[matrix.length][];
        for (int i = 0; i < starsByCol.length;  i++) {
            retval[i] = new int[]{starsByCol[i],i};
        }
        return retval;
    }

    static private boolean allAreCovered(int[] coveredCols) {
        for (int covered : coveredCols) {
            if (0 == covered) return false;
        }
        return true;
    }


    /**
     * the first step of the hungarian algorithm
     * is to find the smallest element in each row
     * and subtract it's values from all elements
     * in that row
     *
     * @return the next step to perform
     */
    static private void reduceMatrix(double[][] matrix) {

        for (int i = 0; i < matrix.length; i++) {

            // find the min value in the row
            double minValInRow = Double.MAX_VALUE;
            for (int j = 0; j < matrix[i].length; j++) {
                if (minValInRow > matrix[i][j]) {
                    minValInRow = matrix[i][j];
                }
            }

            // subtract it from all values in the row
            for (int j = 0; j < matrix[i].length; j++) {
                matrix[i][j] -= minValInRow;
            }
        }

        for (int i = 0; i < matrix[0].length; i++) {
            double minValInCol = Double.MAX_VALUE;
            for (int j = 0; j < matrix.length; j++) {
                if (minValInCol > matrix[j][i]) {
                    minValInCol = matrix[j][i];
                }
            }

            for (int j = 0; j < matrix.length; j++) {
                matrix[j][i] -= minValInCol;
            }

        }

    }

    /**
     * init starred zeroes
     *
     * for each column find the first zero
     * if there is no other starred zero in that row
     * then star the zero, cover the column and row and
     * go onto the next column
     *
     * @param costMatrix
     * @param starredZeroes
     * @param coveredRows
     * @param coveredCols
     * @return the next step to perform
     */
    static private void initStars(double costMatrix[][], int[] starsByRow, int[] starsByCol) {


        int [] rowHasStarredZero = new int[costMatrix.length];
        int [] colHasStarredZero = new int[costMatrix[0].length];

        for (int i = 0; i < costMatrix.length; i++) {
            for (int j = 0; j < costMatrix[i].length; j++) {
                if (0 == costMatrix[i][j] && 0 == rowHasStarredZero[i] && 0 == colHasStarredZero[j]) {
                    starsByRow[i] = j;
                    starsByCol[j] = i;
                    rowHasStarredZero[i] = 1;
                    colHasStarredZero[j] = 1;
                    break; // move onto the next row
                }
            }
        }
    }


    /**
     * just marke the columns covered for any coluimn containing a starred zero
     * @param starsByCol
     * @param coveredCols
     */
    static private void coverColumnsOfStarredZeroes(int[] starsByCol, int[] coveredCols) {
        for (int i = 0; i < starsByCol.length; i++) {
            coveredCols[i] = -1 == starsByCol[i] ? 0 : 1;
        }
    }


    /**
     * finds some uncovered zero and primes it
     * @param matrix
     * @param primesByRow
     * @param coveredRows
     * @param coveredCols
     * @return
     */
    static  private int[] primeSomeUncoveredZero(double matrix[][], int[] primesByRow,
                                       int[] coveredRows, int[] coveredCols) {


        // find an uncovered zero and prime it
        for (int i = 0; i < matrix.length; i++) {
            if (1 == coveredRows[i]) continue;
            for (int j = 0; j < matrix[i].length; j++) {
                // if it's a zero and the column is not covered
                if (0 == matrix[i][j] && 0 == coveredCols[j]) {

                    // ok this is an unstarred zero
                    // prime it
                    primesByRow[i] = j;
                    return new int[]{i,j};
                }
            }
        }
        return null;

    }

    /**
     *
     * @param unpairedZeroPrime
     * @param starsByRow
     * @param starsByCol
     * @param primesByRow
     */
    static  private void incrementSetOfStarredZeroes(int[] unpairedZeroPrime, int[] starsByRow, int[] starsByCol, int[] primesByRow) {

        // build the alternating zero sequence (prime, star, prime, star, etc)
        int i, j = unpairedZeroPrime[1];

        Set<int[]> zeroSequence = new LinkedHashSet<int[]>();
        zeroSequence.add(unpairedZeroPrime);
        boolean paired = false;
        do {
            i = starsByCol[j];
            paired = -1 != i && zeroSequence.add(new int[]{i,j});
            if (!paired) break;

            j = primesByRow[i];
            paired = -1 != j && zeroSequence.add(new int[]{ i, j });

        } while (paired);


        // unstar each starred zero of the sequence
        // and star each primed zero of the sequence
        for (int[] zero : zeroSequence) {
            if (starsByCol[zero[1]] == zero[0]) {
                starsByCol[zero[1]] = -1;
                starsByRow[zero[0]] = -1;
            }
            if (primesByRow[zero[0]] == zero[1]) {
                starsByRow[zero[0]] = zero[1];
                starsByCol[zero[1]] = zero[0];
            }
        }

    }


    static  private void makeMoreZeroes(double[][] matrix, int[] coveredRows, int[] coveredCols) {

        // find the minimum uncovered value
        double minUncoveredValue = Double.MAX_VALUE;
        for (int i = 0; i < matrix.length; i++) {
            if (0 == coveredRows[i]) {
                for (int j = 0; j < matrix[i].length; j++) {
                    if (0 == coveredCols[j] && matrix[i][j] < minUncoveredValue) {
                        minUncoveredValue = matrix[i][j];
                    }
                }
            }
        }

        // add the min value to all covered rows
        for (int i = 0; i < coveredRows.length; i++) {
            if (1 == coveredRows[i]) {
                for (int j = 0; j < matrix[i].length; j++) {
                        matrix[i][j] += minUncoveredValue;
                }
            }
        }

        // subtract the min value from all uncovered columns
        for (int i = 0; i < coveredCols.length; i++) {
            if (0 == coveredCols[i]) {
                for (int j = 0; j < matrix.length; j++) {
                    matrix[j][i] -= minUncoveredValue;
                }
            }
        }
    }
	public static double RatedAssign(int [] L,  double [][] Q, int [][] T, int m, int n, double th) {
		double v=0.0f;
		int cnt=0, LL[]=new int [m];
		double Q1[][]= new double [m][n];
		//Check if it can be a square matrix.
		for (int i = 0; i<n; i++) cnt +=L[i];
		if (cnt > m) return 0.0f;//Not enough agents.
		//Adjust Q with the threshold.
		for (int i = 0; i < m; i++)
			for (int j =0; j< n; j++) 
				if (Q[i][j]<=th) Q1[i][j]=-m*m;
				else Q1[i][j]=Q[i][j];
		double CC[]= new double [n]; //		CC[] is the numbers of qualified agents for roles.
		int D[]= new int [n];		//D is a vector for the difference between the number required current agents and the number of actual current agents for each role.   
		for (int j =0; j< n; j++)
		{	for (int i = 0; i < m; i++)			{		if (Q1[i][j]>th)	 CC[j]=CC[j]+1;			}
			D[j]=(int)( L[j]-CC[j]);	
			//	if (D[j]>0) return 0.0f;//One role has not enough agents..
		}
		//Create the index vector.
		int index=0;
		for (int j = 0; j<n; j++)
				for (int k = 0; k<L[j]; k++) LL[index++] =j;
		for (int k = index; k < m; k++)//fill the left columns (roles).
			LL[k]=n;
//		System.out.println("Q1:");		for (int i = 0; i<m; i++)		{ for (int j = 0; j<n; j++)	System.out.print(Q1[i][j]+" ");		System.out.println();		}		System.out.println();

		double [][] M = new double [m][m];
		for (int i = 0; i<m; i++)
		{ 	index =0;
			for (int j = 0; j<n; j++)
			{
				for (int k = 0; k<L[j]; k++)
					M[i][index++]=1-Q1[i][j];
			}
			//for (int k = index; k < m; k++)//fill the left columns (roles).
				//M[i][k]=1;
		}

		
//	for (int i = 0; i < m; i++)			{System.out.print (LL[i]);System.out.print (" ");	}		System.out.println ();
//	for (int i = 0; i < m; i++)		{	for (int j =0; j< m; j++)			{				System.out.print (M[i][j]);	System.out.print (" ");	}System.out.println ();		}
		int [][] N = computeAssignments(M);
		//Obtaing the matrix T.
		for (int i = 0; i<m; i++)
		{ 	for (int j = 0; j<n; j++)
			{					T[i][j]=0;
			}
		}
		
		for (int i = 0; i < N.length; i++)
			if (LL[N[i][1]]< n)
			{
				if (Q1[N[i][0]][LL[N[i][1]]]>0.0f) T[N[i][0]][LL[N[i][1]]]=1;
				else return -1.0f;
				T[N[i][0]][LL[N[i][1]]]=1;
			}
	for (int i = 0; i<m; i++)
		 	for (int j = 0; j<n; j++)
				v += Q[i][j]*T[i][j];
		return v;
	}
}
class MatrixIndex{
	int i, j;
	MatrixIndex(int x, int y){
		i = x;
		j = y;
	}
}

public class MERAGRA {
	private static String filename="Result.txt";
	static DecimalFormat tw = new DecimalFormat("0.00");
	 public static  void printDMatrix (BufferedWriter out, double [][]x, int m, int n){
	  		try 
			{
			for (int i = 0; i < m; i++)
			{	for (int j =0; j< n; j++)
				{
				out.write (tw.format(x[i][j]));		out.write (" ");
				}
			out.write ("\n");
			}
			out.write ("\n");
			}
			catch (IOException e) {System.out.println ("Error in writing into a file!");}	
		}	
	 public static  void printIMatrix (BufferedWriter out, int [][]x, int m, int n){
			DecimalFormat tw = new DecimalFormat("0");
	  		try 
			{
			for (int i = 0; i < m; i++)
			{	for (int j =0; j< n; j++)
				{
				out.write (tw.format(x[i][j]));		out.write (" ");
				}
			out.write ("\n");
			}
			out.write ("\n");
			}
			catch (IOException e) {System.out.println ("Error in writing into a file!");}	
		}	
	 public static  void printDMatrix (double [][]x, int m, int n){
			DecimalFormat tw = new DecimalFormat("0.00");
			for (int i = 0; i < m; i++)
			{	for (int j =0; j< n; j++)
				{
				System.out.print (tw.format(x[i][j]));		System.out.print (" ");
				}
			System.out.println ();
			}
			System.out.println ();
		}	
	 public static  void printIMatrix (int [][]x, int m, int n){
			DecimalFormat tw = new DecimalFormat("0");
			for (int i = 0; i < m; i++)
			{	for (int j =0; j< n; j++)
				{
				System.out.print (tw.format(x[i][j]));		System.out.print (" ");
				}
			System.out.println ();
			}
			System.out.println ();
		}
	public static int maxIndex(double [][]V, int j){
		double max = 0;
		int index=0;
		for (int i = 0; i < V.length; i++){
			if (V[i][j]>max){max=V[i][j]; index=i;}
		}
		return index;
	}
	public static int sumRowTo(int []V, int k){
		int sum = 0;
		for (int j = 0; j < k; j++)	sum+=V[j];
		return sum;
	}
	public static int NoteRoleInRow(int [] T){
		for (int j = 0; j < T.length; j++) if (T[j]==1) return j;
		return 0;
	}
	public static boolean Satisfied(int [][]T, int []L, MatrixIndex md, int m, int n){
		boolean res=true;
		int sum =0; 
		int role=0;
		for (int j = 0; j < n; j++) sum+=T[md.i][j];
		if (sum>=2) return false;
		if (sum==1) {
			for (int j = 0; j < n; j++) if (T[md.i][j]==1) role = j;
			for (int i = 0; i < m; i++){
				if ((T[i][role]==1)&&(T[i][md.j]==1)) return false;
			}
		}
		sum =0;
		for (int i = 0; i < m; i++) sum+=T[i][md.j];
		if (sum==(L[md.j]+1)) return false;
		return res;
	}	 
		public static boolean GMEO_GRA(int m, int n, double [][]Q, int []L, int [][]QD,int na){
			boolean res = true;
			double Q1[][]=new double [m][n];
			int nt=na+n;
			if (nt>m) return false;// No solution!
			int [][]T=new int [m][n];
			RatedMunkres.RatedAssign(L, Q, T, m, n, 0);
			for (int i = 0; i<m;i++) for (int j =0; j< n; j++) {
				if (T[i][j]!=1) Q1[i][j]=Q[i][j];
				else {Q1[i][j]=-m*n;}
			}
			Set<Integer> LeftRSet = new HashSet<Integer>();
			for (int j =0; j< n; j++) LeftRSet.add(j);
			MatrixIndex md=new MatrixIndex(0,0);
			int []TagA=new int [m];
			for (int i = 0; i<m;i++) TagA[i]=0; 
			double kmax=0;
			for (int i = 0; i<m;i++){ 
				int ss=0;
				for (int j =0; j< n; j++) {if (T[i][j]==0) ss++;}
				if (ss==n) { TagA[i]=1;
					for (int j =0; j< n; j++)
						if (Q1[i][j]> kmax) {md.i=i; md.j=j;kmax=Q1[i][j];}
				}
			}			
			T[md.i][md.j]=1;
			Q1[md.i][md.j]=-m*n;
			LeftRSet.remove(md.j);	
			do {
				kmax=0;
				for (int i = 0; i<m;i++){ 
					for (int j:LeftRSet) 
						if (Q1[i][j]> kmax) {md.i=i; md.j=j;kmax=Q1[i][j];}
				}
				if (Satisfied(T, L, md, m, n)){
					T[md.i][md.j]=1;
					Q1[md.i][md.j]=-m*n;			
					LeftRSet.remove(md.j);
				}
				else {
					Q1[md.i][md.j]=-m*n;
				}
			} while (!LeftRSet.isEmpty());
			for (int i = 0; i<m;i++) for (int j =0; j< n; j++) QD[i][j]=T[i][j];
			return res;
		}
	public static void main(String[] args)
	{	
		Random generator = new Random();
		DecimalFormat df = new DecimalFormat("0.00");
		
		int m = 13, n = 4;
		int L[]= {1, 2, 4, 2};
		int []VC = {3,5,2,4};
		double tao = 0.85; 
		double [][]Q={
				{0.18,0.82,0.29,0.01},
				{0.35,0.80,0.58,0.35},
				{0.84,0.85,0.86,0.36},
				{0.96,0.51,0.45,0.64},
				{0.22,0.33,0.68,0.33},
				{0.96,0.50,0.10,0.73},
				{0.25,0.18,0.23,0.39},
				{0.56,0.35,0.80,0.62},
				{0.49,0.09,0.33,0.58},
				{0.38,0.54,0.72,0.20},
				{0.91,0.31,0.34,0.15},
				{0.85,0.34,0.43,0.18},
				{0.44,0.06,0.66,0.37}};
		int[][] QD = new int[m][n];
		int[][] TD = new int[m][n];
		long t1 = System.nanoTime();
		int na = 0;
		for (int j =0; j<n; j++) {na+=L[j];};	
		boolean res=GMEO_GRA(m, n, Q, L, QD, na);
		long t2 = System.nanoTime();
		long diff = (t2-t1)/1000000;
		printIMatrix(QD, m, n);
		double v = 0;
		for(int r=0; r<m; r++)		{
			for(int c=0; c<n; c++)			{
				if ((QD[r][c] ==1)&&(Q[r][c]<tao)) TD[r][c] =1;
				else TD[r][c] =0;
				if (QD[r][c] ==1) v+=Q[r][c];
			}
		}
		int Cost =0;
		int NL[]=new int [n];
			for(int c=0; c<n; c++)			{
				NL[c]=0;
				for(int r=0; r<m; r++)	NL[c]+= TD[r][c];
				Cost+=NL[c]*VC[c];
		}
			System.out.println("T= ");
			printIMatrix(QD, m, n);
			System.out.println("TD = ");
			printIMatrix(TD, m, n);
			System.out.println("time= "+diff+"ms");
			System.out.println("v= "+df.format(v));
			System.out.println("cost= "+Cost);		
	}
}