//
// assign - computes a maximum weighted bipartite b-matching between
//          persons and tasks
//
// Copyright (C) 2003-2020  Joachim Reichel <joachim.reichel@posteo.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 <cassert>
#include <climits>
#include <cstring>
#include <fstream>
#include <iostream>
#include <sstream>
#include <vector>

#ifdef USE_CPLEX
#include "ilcplex/ilocplex.h"
#else // USE_CPLEX
#include "coin/CbcModel.hpp"
#include "coin/CoinModel.hpp"
#include "coin/OsiClpSolverInterface.hpp"
#endif // USE_CPLEX


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

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

std::string get_row_name_person( int p)
{
    std::ostringstream s;
    s << "row_p" << p;
    return s.str();
}

std::string get_row_name_task( int t)
{
    std::ostringstream s;
    s << "row_t" << t;
    return s.str();
}

#ifdef USE_CPLEX
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 "";
}
#endif // USE_CPLEX

int main( int argc, char* argv[]) // NOLINT bugprone-exception-escape
{
    // 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]);

    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+8+t]);

    // 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();

    // 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( const auto& x: wish)
        for( int w = 0; w < nr_of_wishes; ++w)
            wishes_per_task[ x[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;

#ifdef USE_CPLEX

    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_column_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, get_row_name_person( p).c_str());
        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, get_row_name_task( t).c_str());
        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();

#else // USE_CPLEX

    CoinModel model;

    // generate columns

    for( int p = 0; p < nr_of_persons; ++p)
        for( int t = 0; t < nr_of_tasks; ++t) {
            int i = p*nr_of_tasks + t;
            model.setColumnName( i, get_column_name( p, t).c_str());
            model.setColumnBounds( i, 0, 1);
            model.setInteger( i);
        }

    // generate objective

    for( int p = 0; p < nr_of_persons; ++p)
        for( int t = 0; t < nr_of_tasks; ++t)
            model.setObjective( 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)
            model.setObjective( p*nr_of_tasks + wish[p][w], coeff[w]);

    // generate rows

    for( int p = 0; p < nr_of_persons; ++p) {
        model.setRowName( p, get_row_name_person( p).c_str());
        model.setRowBounds( p, tpp_min, tpp_max);
        for( int t = 0; t < nr_of_tasks; ++t)
            model.setElement( p, p*nr_of_tasks + t, 1);
    }

    for( int t = 0; t < nr_of_tasks; ++t) {
        int i = nr_of_persons + t;
        model.setRowName( i, get_row_name_task( t).c_str());
        model.setRowBounds( i, stpt[t] * ppt_min, stpt[t] * ppt_max);
        for( int p = 0; p < nr_of_persons; ++p)
            model.setElement( i, p*nr_of_tasks + t, 1);
    }

    // workaround for a bug

    model.setRowBounds( 0, tpp_min, tpp_max);
    assert( model.getRowLower( 0) == tpp_min);
    assert( model.getRowUpper( 0) == tpp_max);

    // solve model

    OsiClpSolverInterface solver;
    solver.loadFromCoinModel( model);

    CbcModel cbc_model( solver);
    cbc_model.setObjSense( -1);
    cbc_model.solver()->setHintParam( OsiDoReducePrint, true, OsiHintTry);
    cbc_model.messageHandler()->setLogLevel( 0);
    cbc_model.branchAndBound();

    if( cbc_model.status()) {
        std::cout << "No optimal solution found" << std::endl;
        return 2;
    }

    std::vector<double> x( nr_of_columns);
    const double* tmp = cbc_model.solver()->getColSolution();
    for( int i = 0; i < nr_of_columns; ++i)
        x[i] = tmp[i];

    double obj = cbc_model.getObjValue();

#endif // USE_CPLEX

    // 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::vector<int>> p2t( nr_of_persons);
    std::vector<std::vector<int>> t2p( nr_of_tasks);

    for( int i = 0; i < nr_of_columns; ++i) {

        if( x[i] == 0)
            continue;

        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( const auto& x: p2t[p])
            std::cout << x+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( const auto& x: t2p[t])
             std::cout << x+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;
}