/*
		$Id: assign.cpp 3343 2008-02-19 15:31:56Z joachim $

		assign - computes a maximum weighted bipartite b-matching between
		         persons and tasks


		Copyright (C) 2003, 2008  Joachim Reichel <joachim.reichel@gmx.de>

		Redistribution and use in source and binary forms, with or without
		modification, are permitted provided that the following conditions
		are met:

		1. Redistributions of source code must retain the above copyright
		   notice, this list of conditions and the following disclaimer.
		2. Redistributions in binary form must reproduce the above copyright
		   notice, this list of conditions and the following disclaimer in the
		   documentation and/or other materials provided with the distribution.
		3. Neither the name of the University nor the names of its contributors
		   may be used to endorse or promote products derived from this software
		   without specific prior written permission.

		THIS SOFTWARE IS PROVIDED BY THE REGENTS AND CONTRIBUTORS ``AS IS'' AND
		ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
		IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
		ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR CONTRIBUTORS BE LIABLE
		FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
		DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
		OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
		HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
		LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
		OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
		SUCH DAMAGE.
*/

/*
		The program computes a maximum weighted bipartite b-matching between
		persons and tasks. Possible use cases are
		(1) assign students (persons) to tutorial groups (tasks)
		(2) assign students (persons) to seminar topics (tasks)
		(3) assign advisors (persons) to seminar topics (tasks)

		The implementation is able to solve the general weighted bipartite
		b-matching problem. However, the problem specification places restrictions
		on the weights. Furthermore, only uniform lower and upper bounds for the
		components of the b-vector can be specified.

		The persons are identified by nodes from 1 to p, the tasks by nodes
		from 1 to t. Lower and upper bounds on the number of tasks per person are
		given by tpp_min and tpp_max. Similarly, lower and upper bounds on the
		number of persons per task are given by ppt_min and ppt_max. In the above
		use cases, these bounds should be set as follows:
		(1) tpp_min = tpp_max = 1, ppt_min and ppt_max as desired
		(2) tpp_min = tpp_max = 1, ppt_min = ppt_max = 1
		(3) ppt_min = ppt_max = 1, tpp_min and tpp_max as desired
		In case (2) tpp_min should be set to 0 if there are more students than
		seminar topics. If there are less students than topics, set ppt_min to 0.

		The graph is the complete bipartite graph with weights as follows. Each
		person expresses w preferences for certain tasks. The edge connecting
		a person node with the task node corresponding to its i-th preference
		receives weight c_i. All other edges receive weight c_0 (which is usually
		much smaller than c_w).

		The values n_i can be used to group several indistinguishable tasks
		(called subtasks). For example, two tutorial groups at the same time slot
		are represented by *one* task with two subtasks (assuming that students
		do not select a particular tutorial group but a time slot.). The value n_i
		indicates the number of subtasks of task i and defaults to 1. (In the
		implemention, the lower and upper bounds for task node i are determined by
		the product of ppt_min and n_i and ppt_max and n_i, respectively.)

		The file format is as follows: There is one line per person, the number
		p of persons is implicitly given by the number of lines. On each line,
		there are w integers from 1 to t indicating the preferences of that
		person. The i-th entry corresponds to the i-th preference.

		Note 1: Often there are several optimal solutions. However, only one such
		solution is computed. Keep in mind that slight variations in the
		objective coefficients can have a large influence.

		Note 2: If a person expresses less than w preferences, simply repeat
		the lowest preferences as often as needed (assuming that the c_i,
		1 <= i <= w, are decreasing).

		Warning: There is no error checking for the input.
*/



#include <iostream>
#include <fstream>
#include <vector>
#include <list>

#include "ilcplex/ilocplex.h"



void usage ()
{
	std::cout << "Usage: assign <file> <t> <w> <tpp_min> <tpp_max> <ppt_min ppt_max> \\" << std::endl;
	std::cout << "                <c_1> ... <c_w> <c_0> [<n_1> <n_2> ... <n_t>]" << std::endl;
}


std::string get_name (int p, int t)
{
	std::ostringstream S;
	S << "x_p" << p << "_t" << t;
	return S.str();
}


std::string cplex_status (IloCplex& cplex)
{
	switch (cplex.getStatus())
		{
			case IloAlgorithm::Unknown:               return "unknown status";
			case IloAlgorithm::Feasible:              return "feasible";
			case IloAlgorithm::Optimal:               return "optimal";
			case IloAlgorithm::Infeasible:            return "infeasible";
			case IloAlgorithm::Unbounded:             return "unbounded";
			case IloAlgorithm::InfeasibleOrUnbounded: return "infeasible or unbounded";
			case IloAlgorithm::Error:                 return "error";
			default:                                  assert (false);
		}
	return "";
}


int main (int argc, char* argv[])
{
	// read command line arguments

	if (argc < 1+6)
		{
			usage ();
			return 1;
		}

	std::string filename = argv[1];
	int nr_of_tasks = atoi (argv[2]);
	int nr_of_wishes = atoi (argv[3]);

	if ((argc != 1+nr_of_wishes+8) && (argc != 1+nr_of_wishes+8+nr_of_tasks))
		{
			usage ();
			return 1;
		}

	int tpp_min = atoi (argv[4]);
	int tpp_max = atoi (argv[5]);
	int ppt_min = atoi (argv[6]);
	int ppt_max = atoi (argv[7]);

	std::vector<double> coeff (nr_of_wishes+1);
	for (int w=0; w<=nr_of_wishes; w++)
		coeff[w] = atof (argv[w+8]);

	// read data file

	std::vector<std::vector<int> > wish;
	int nr_of_persons = 0;
	std::ifstream file (filename.c_str());
	while (!file.eof())
		{
			std::vector<int> tmp (nr_of_wishes);
			for (int w=0; w<nr_of_wishes; w++)
				{
					file >> tmp[w];
					if (file.eof())
						break;
					tmp[w]--;
				}
			if (file.eof())
				break;
			wish.push_back (tmp);
			nr_of_persons++;
		}
	file.close();

	std::vector<int> stpt (nr_of_tasks, 1);
	if (argc > 1+nr_of_wishes+8)
		for (int t=0; t<nr_of_tasks; t++)
			stpt[t] = atoi (argv[1+nr_of_wishes+6+t]);

	// print input

	std::cout << "filename: " << filename << std::endl;
	std::cout << "# persons: " << nr_of_persons << std::endl;
	std::cout << "# tasks: " << nr_of_tasks << std::endl;
	std::cout << "# wishes: " << nr_of_wishes << std::endl;
	std::cout << "min. #tasks/person: " << tpp_min << std::endl;
	std::cout << "max. #tasks/person: " << tpp_max << std::endl;
	std::cout << "min. #persons/task: " << ppt_min << std::endl;
	std::cout << "max. #persons/task: " << ppt_max << std::endl;
	std::cout << std::endl;

	for (int w=0; w<nr_of_wishes; w++)
		std::cout << "obj. coeff. for " << w+1 << ". wish fulfilled: " << coeff[w] << std::endl;
	std::cout << "obj. coeff. for no wish fulfilled: " << coeff[nr_of_wishes] << std::endl;
	std::cout << std::endl;

	for (int t=0; t<nr_of_tasks; t++)
		std::cout << "# subtasks for task " << t+1 << ": " << stpt[t] << std::endl;
	std::cout << std::endl;

	// some statistics

	std::vector<std::vector<int> > wishes_per_task (nr_of_tasks);
	for (int t=0; t<nr_of_tasks; t++)
		wishes_per_task[t].resize (nr_of_wishes+1, 0);
	for (std::vector<std::vector<int> >::iterator i = wish.begin(); i != wish.end(); i++)
		for (int w=0; w<nr_of_wishes; w++)
			wishes_per_task[ (*i)[w] ][w]++;

	 for (int t=0; t<nr_of_tasks; t++)
		{
			std::cout << "# wishes for task " << t+1 << ": ";
			for (int w=0; w<nr_of_wishes; w++)
				std::cout << wishes_per_task[t][w] << " ";
			std::cout << std::endl;
		}
	std::cout << std::endl;

	// setup ILP

	int nr_of_rows = nr_of_persons + nr_of_tasks;
	int nr_of_columns = nr_of_persons * nr_of_tasks;
	int nr_of_nonzeros = 2 * nr_of_persons * nr_of_tasks;

	std::cout << "nr_of_rows     : " << nr_of_rows << std::endl;
	std::cout << "nr_of_columns  : " << nr_of_columns << std::endl;
	std::cout << "nr_of_nonzeros : " << nr_of_nonzeros << std::endl;
	std::cout << std::endl;

	IloEnv env;
	IloCplex cplex (env);
	cplex.setOut (env.getNullStream());

	IloModel model (env);

	// generate columns

	IloIntVarArray vars (env);
	for (int p=0; p<nr_of_persons; p++)
		for (int t=0; t<nr_of_tasks; t++)
			{
				IloIntVar var (env, 0, 1, get_name (p, t).c_str());
				vars.add (var);
			}
	model.add (vars);

	// generate objective

	IloObjective objective = IloMaximize (env);
	for (int p=0; p<nr_of_persons; p++)
		for (int t=0; t<nr_of_tasks; t++)
			objective.setLinearCoef (vars[p*nr_of_tasks + t], coeff[nr_of_wishes]);
	for (int p=0; p<nr_of_persons; p++)
		for (int w=0; w<nr_of_wishes; w++)
			objective.setLinearCoef (vars[p*nr_of_tasks + wish[p][w]], coeff[w]);
	model.add (objective);

	// generate rows

	for (int p=0; p<nr_of_persons; p++)
		{
			IloRange range (env, tpp_min, tpp_max);
			for (int t=0; t<nr_of_tasks; t++)
				range.setLinearCoef (vars[p*nr_of_tasks + t], 1);
			model.add (range);
		}

	for (int t=0; t<nr_of_tasks; t++)
		{
			IloRange range (env, stpt[t] * ppt_min, stpt[t] * ppt_max);
			for (int p=0; p<nr_of_persons; p++)
				range.setLinearCoef (vars[p*nr_of_tasks + t], 1);
			model.add (range);
		}

	// solve model

	cplex.extract (model);
	cplex.solve();

	if (cplex.getStatus() != IloAlgorithm::Optimal)
		{
			std::cout << "No optimal solution found, reason: " << cplex_status (cplex) << std::endl;
			return 2;
		}

	IloNumArray tmp (env);
	cplex.getValues (tmp, vars);

	std::vector<double> x (nr_of_columns);
	for (int i=0; i<nr_of_columns; i++)
		x[i] = tmp[i];
	tmp.end();

	double obj = cplex.getObjValue();

	// destroy model

	IloExtractableArray remove (env);
	for (IloModel::Iterator it (model); it.ok(); ++it)
		remove.add (*it);
	remove.add (model);
	remove.endElements();
	remove.end();
	vars.end();

	// output

	std::vector<int> tpp_size (nr_of_persons, 0);
	std::vector<int> ppt_size (nr_of_tasks, 0);
	std::vector<int> wishes_fulfilled (nr_of_wishes+1, 0);

	std::vector<std::list<int> > p2t (nr_of_persons);
	std::vector<std::list<int> > t2p (nr_of_tasks);

	for (int i=0; i<nr_of_columns; i++)
		if (x[i] != 0)
			{
				int person = i/nr_of_tasks;
				int task = i%nr_of_tasks;

				tpp_size[person]++;
				ppt_size[task]++;
				p2t[person].push_back (task);
				t2p[task].push_back (person);

				bool fulfilled = false;
				for (int w=0; w<nr_of_wishes; w++)
					if (wish[person][w] == task)
						{
							fulfilled = true;
							wishes_fulfilled[w]++;
							break;
						}
				if (!fulfilled)
					wishes_fulfilled[nr_of_wishes]++;
			}

	for (int p=0; p<nr_of_persons; p++)
		{
			std::cout << "person " << p+1 << " assigned to task(s) ";
			for (std::list<int>::iterator it = p2t[p].begin(); it != p2t[p].end(); it++)
				std::cout << (*it)+1 << " ";
			 std::cout << std::endl;
		}
	std::cout << std::endl;

	for (int t=0; t<nr_of_tasks; t++)
		{
			std::cout << "task " << t+1 << " assigned to person(s) ";
			for (std::list<int>::iterator it = t2p[t].begin(); it != t2p[t].end(); it++)
				std::cout << (*it)+1 << " ";
			 std::cout << std::endl;
		}
	std::cout << std::endl;

	std::cout << "objective value: " << obj << std::endl;
	std::cout << std::endl;

	for (int p=0; p<nr_of_persons; p++)
		std::cout << "tasks for person " << p+1 << ": " << tpp_size[p] << std::endl;
	std::cout << std::endl;

	for (int t=0; t<nr_of_tasks; t++)
		std::cout << "persons for task " << t+1 << ": " << ppt_size[t] << std::endl;
	std::cout << std::endl;

	for (int w=0; w<nr_of_wishes; w++)
		std::cout << "# " << w+1 << ". wish fulfilled: " << wishes_fulfilled[w] << std::endl;
	std::cout << "# no wish fulfilled: " << wishes_fulfilled[nr_of_wishes] << std::endl;
	std::cout << std::endl;

	return 0;
}