#include <ilcplex/ilocplex.h>
#include <vector>

ILOSTLBEGIN

typedef IloArray <IloNumVarArray> IloNumVarArray2;

int sum(int num)
{
	int totalnum = 0;
	
	for(int x1=1; x1<=num; x1++)
	{
		totalnum += x1;
	}
	
	return totalnum;
}

int assign_index(int i, int j)
{

	int k;
	
	if(i==0 || j==0)
	{
		return -999;
	}
	
	else
	{
		k = sum(i) - (i-j);
	}


	return k;
}



int main()
{
	

	try {

		ofstream wf;

		const int N = 64;		//	The number of users
		const int C = N;		//	Network capacity
		
		//double lambda[] = {0.01, 0.02, 0.04, 0.1, 0.15, 0.2, 0.34, 0.5, 1, 1.34, 2, 4, 10, 16.67, 20, 33.34, 50, 100}; 
		//double lambda[] = {50, 100}; 
		double mu = 1.0;
		vector<double> lambda;	//	Arrival rates
		lambda.push_back(0.01);
		lambda.push_back(0.02);
		lambda.push_back(0.04);
		lambda.push_back(0.1);
		lambda.push_back(0.15);
		lambda.push_back(0.2);
		lambda.push_back(0.34);
		lambda.push_back(0.5);
		lambda.push_back(1);
		lambda.push_back(1.34);
		lambda.push_back(2);
		lambda.push_back(4);
		lambda.push_back(10);
		lambda.push_back(16.67);
		lambda.push_back(20);
		lambda.push_back(33.34);
		lambda.push_back(50);
		lambda.push_back(100);
		

		for(size_t AR=0; AR<lambda.size(); AR++){
			//cout << "Lambda: " << lambda[AR] << endl;

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

			//wf.open("blocking_newmodel_ILP_N64_C32_test2.txt", ios::out | ios::app);
			wf.open("Blocking_newmodel_ILP_N64_C32_test2.txt", ios::out | ios::app);

			/***** P *****/
			IloNumVarArray2 P(env, N+1);
			for(int i1=0; i1 <= N; i1++)
			{
				P[i1] = IloNumVarArray(env, N+1);
		
				for(int j1=0; j1 <= N; j1++)
				{
					P[i1][j1] = IloNumVar(env, 0.0, 1.0, IloNumVar::Float);
					//cout << " i: " << i << " j:" << j << endl; 
				}
			}

			



			/*****	Constraints	*****/
			for(int i=0; i<=N; i++)
			{
				for(int j=0; j<=N; j++)
				{
					//cout << "i: " << i << " j:" << j << endl;
					if(i==0 && j==0)	//	P00 
					{
						model.add( (i+1)*mu * P[i+1][j+1] == (N-j)*lambda[AR]*(N-i)/N * P[i][j] );
						//model.add( (j+1)*mu * P[i+1][j+1] == (N-i)*lambda[AR]*(N-j)/N * P[i][j] );
						//wf << model.add( (i+1)*mu * P[i+1][j+1] == (N-j)*lambda[AR]*(N-i)/N * P[i][j] ) << endl;
					}

					else if(i==1 && j==1)	//	P11
					{
						model.add( (N-j+1)*lambda[AR]*(N-i+1)/N * P[i-1][j-1] + i*mu*(j-i+1)/j * P[i][j+1] + (i+1)*mu*i/j * P[i+1][j+1] 
								   == ( (i*mu*i/j) + (N-j)*lambda[AR] ) * P[i][j]  );
						//model.add( (N-i+1)*lambda[AR]*(N-j+1)/N * P[i-1][j-1] + j*mu*(i-j+1)/i * P[i+1][j] + (j+1)*mu*j/i * P[i+1][j+1] 
						//		   == ( (j*mu*j/i) + (N-i)*lambda[AR] ) * P[i][j]  );
					}

					else if(i==1 && j>=2 && j<N)	//	upper boundary condition (P21, P31, P41...)
					{
						model.add( (N-j+1)*lambda[AR]*i/N * P[i][j-1] + i*mu*(j-i+1)/j * P[i][j+1] + (i+1)*mu*i/j * P[i+1][j+1]
									== ( i*mu*(j-i)/(j-1) + (N-j)*lambda[AR] ) * P[i][j]);
						//model.add( (N-i+1)*lambda[AR]*j/N * P[i-1][j] + j*mu*(i-j+1)/i * P[i+1][j] + (j+1)*mu*j/i * P[i+1][j+1]
						//			== ( j*mu*(i-j)/(i-1) + (N-i)*lambda[AR] ) * P[i][j]);
					}

					else if(i==1 && j==N)		//	The end of upper boundary (P_N1)
					{
						model.add( (N-j+1)*lambda[AR]*i/N * P[i][j-1] == i*mu*(j-i)/(j-1) * P[i][j] );
						//model.add( (N-i+1)*lambda[AR]*j/N * P[i-1][j] == j*mu*(i-j)/(i-1) * P[i][j] );
					}

					else if(i==j && j>=2 && j<N)	//	diagonal condition (P22, P33, P44...)
					{
						if(i<C)
						{
							model.add( (N-j+1)*lambda[AR]*(N-i+1)/N * P[i-1][j-1] + i*mu*(j-i+1)/j * P[i][j+1] + (i+1)*mu*i/j * P[i+1][j+1]
										== ( i*mu*(i-1)/(j-1) + (N-j)*lambda[AR] ) * P[i][j] );
							//model.add( (N-i+1)*lambda[AR]*(N-j+1)/N * P[i-1][j-1] + j*mu*(i-j+1)/i * P[i+1][j] + (j+1)*mu*j/i * P[i+1][j+1]
							//			== ( j*mu*(j-1)/(i-1) + (N-i)*lambda[AR] ) * P[i][j] );
						}

						else	//	i >= C
						{
							model.add( (N-j+1)*lambda[AR]*(N-i+1)/N * P[i-1][j-1] + C*mu*(j-i+1)/j * P[i][j+1] + C*mu*i/j * P[i+1][j+1]
								   == ( C*mu*(i-1)/(j-1) + (N-j)*lambda[AR] ) * P[i][j] );
							//model.add( (N-i+1)*lambda[AR]*(N-j+1)/N * P[i-1][j-1] + C*mu*(i-j+1)/i * P[i+1][j] + C*mu*j/i * P[i+1][j+1]
							//	   == ( C*mu*(j-1)/(i-1) + (N-i)*lambda[AR] ) * P[i][j] );
						}		
					}

					else if(j==N && i>=2 && i<N)	//	last vertical line (P_N2, P_N3, ..., P_NN-1)
					{
						if(i<C)
						{
							model.add( (N-j+1)*lambda[AR]*(N-i+1)/N * P[i-1][j-1] + (N-j+1)*lambda[AR]*i/N * P[i][j-1] == i*mu*P[i][j]);
							//model.add( (N-i+1)*lambda[AR]*(N-j+1)/N * P[i-1][j-1] + (N-i+1)*lambda[AR]*j/N * P[i-1][j] == j*mu*P[i][j]);
						}

						else	//	i >= C
						{
							model.add( (N-j+1)*lambda[AR]*(N-i+1)/N * P[i-1][j-1] + (N-j+1)*lambda[AR]*i/N * P[i][j-1] == C*mu*P[i][j]);
							//model.add( (N-i+1)*lambda[AR]*(N-j+1)/N * P[i-1][j-1] + (N-i+1)*lambda[AR]*j/N * P[i-1][j] == C*mu*P[i][j]);
						}
					}

					else if(i==N && j==N)
					{
						if(i<C)
						{
							model.add( (N-j+1)*lambda[AR]*(N-i+1)/N * P[i-1][j-1] == i*mu*(i-1)/(j-1) * P[i][j]);
							//model.add( (N-i+1)*lambda[AR]*(N-j+1)/N * P[i-1][j-1] == j*mu*(j-1)/(i-1) * P[i][j]);
						}
						else	//	i >= C
						{
							model.add( (N-j+1)*lambda[AR]*(N-i+1)/N * P[i-1][j-1] == C*mu*(i-1)/(j-1) * P[i][j]);
							//model.add( (N-i+1)*lambda[AR]*(N-j+1)/N * P[i-1][j-1] == C*mu*(j-1)/(i-1) * P[i][j]);
						}
					}

					else if((j==0 && i!=0) || (j!=0 && i==0) || (j<i) )
					{
						model.add(P[i][j] == 0);
					}

					else
					{
						if(i<C)
						{
							model.add( (N-j+1)*lambda[AR]*(N-i+1)/N * P[i-1][j-1] + (N-j+1)*lambda[AR]*i/N * P[i][j-1]
								   + i*mu*(j-i+1)/j * P[i][j+1] + (i+1)*mu*i/j * P[i+1][j+1] == ( i*mu + (N-j)*lambda[AR] ) * P[i][j]);
							//model.add( (N-i+1)*lambda[AR]*(N-j+1)/N * P[i-1][j-1] + (N-i+1)*lambda[AR]*j/N * P[i-1][j]
							//	   + j*mu*(i-j+1)/i * P[i+1][j] + (j+1)*mu*j/i * P[i+1][j+1] == ( j*mu + (N-i)*lambda[AR] ) * P[i][j]);
						}
						
						else	//	i >= C
						{
							model.add( (N-j+1)*lambda[AR]*(N-i+1)/N * P[i-1][j-1] + (N-j+1)*lambda[AR]*i/N * P[i][j-1]
								   + C*mu*(j-i+1)/j * P[i][j+1] + C*mu*i/j * P[i+1][j+1] == ( C*mu + (N-j)*lambda[AR] ) * P[i][j]);
							//model.add( (N-i+1)*lambda[AR]*(N-j+1)/N * P[i-1][j-1] + (N-i+1)*lambda[AR]*j/N * P[i-1][j]
							//	   + C*mu*(i-j+1)/i * P[i+1][j] + C*mu*j/i * P[i+1][j+1] == ( C*mu + (N-i)*lambda[AR] ) * P[i][j]);
						}
						
					}
				}//j
			}//i

			/**	Sum of probability == 1	**/
			IloExpr Psum(env);
			for(int i2=0; i2<=N; i2++)
			{
				for(int j2=0; j2<=N; j2++)
				{
					Psum += P[i2][j2];		
				}//j2
			}//i2
			model.add(Psum == 1);
			//wf << model.add(Psum == 1) << endl;
			Psum.end();

			/**	Range of probability O and 1	**/
			for(int i3=0; i3<=N; i3++)
			{
				for(int j3=0; j3<=N; j3++)
				{
					model.add(P[i3][j3] >= 0);
					//wf << model.add(P[i3][j3] >= 0) << endl;
				}//j3
			}//i3

			for(int i4=0; i4<=N; i4++)
			{
				for(int j4=0; j4<=N; j4++)
				{
					model.add(P[i4][j4] <= 1.0);
					//wf << model.add(P[i4][j4] <= 1.0) << endl;
				}//j4
			}//i4

			/** Throughput **/
			IloExpr Thsum(env);
			for(int i5=0; i5<=N; i5++)
			{
				for(int j5=0; j5<=N; j5++)
				{
					if(i5<C)
					{
						Thsum += i5 * mu * P[i5][j5] / N;
					}
					else
					{
						Thsum += C * mu * P[i5][j5] / N;
					}	
				}//j5
			}//i5
			model.add(IloMaximize(env, Thsum));
			//Thsum.end();

			cplex.solve();

			//wf << "[SOLUTION_STATUS] : " << cplex.getStatus() << endl;
			//wf.precision(4);
			wf.width(7);
			wf.setf(ios::left);
			wf << lambda[AR];
			
			wf.fixed;
			wf.precision(4);
			wf << cplex.getObjValue() << endl;


			wf << endl;
			for(int i=0; i<=N; i++)
			{
				for(int j=0; j<=N; j++)
				{
					wf << setw(13) << cplex.getValue(P[i][j]) << " " ;
				}
				wf << endl;
			}
			wf << endl;

			
			wf.close();
			env.end();

		}//AR

		//cout << "P[0][0]: " << scientific << cplex.getValue(P[0][0]) << endl;
		

		} catch (IloException& e) {
		cerr << "Concert Exception: " << e << endl;
	} 


}//main

