/*                    Linear Assignment Problem
 *
 * Given a set of n agents and n tasks, assign tasks to agent such that
 * every agent has a unique task --- and, of course, every task has
 * someone performing it.  Find the optimal such assignment --- for our
 * purposes here, the one with the largest total benefit.
 *
 * There is a global matrix benefit[MAX_SIZE][MAX_SIZE] which contains
 * for each value of row and col the benefit of assigning agent row to
 * task col.  Because of how C behaves, this has been set to a maximum
 * size, given in a #define.  Also at the global level is a vector to
 * receive the solution and a variable to hold the value of the solution
 * --- this is updated as the backtracking discovered permutations with
 * a better value that the presently stored permutation.
 *
 * The initial solution is generated as one of the two trivial solutions
 * --- either the diagonal solution (0, 1, 2, ...) or the antidiagonal
 * solution (...,2, 1, 0), whichever has the higher value --- as the sum
 * over k of benefit[k][solution[k]].  As a known solution, this
 * provides a lower limit when considering potential permutations.
 *
 * One can also quickly obtain an upper limit on a partial permutation.
 * While we are positioning values as [index], the known value for the
 * permutation will be the sum up to index of benefit[k][solution[k]].
 * Rows from [index+1] to [size-1] have not been assigned yet, and
 * correspondingly columns from perm[index+1] to perm[size-1] have not
 * been assigned.  One can obtain the column maximum value for each of
 * those columns and add them together.  The complete permutation cannot
 * possibly have a value greater than the total of its fixed part and
 * these column maxima, though it well might have one less.  If this
 * number is less than the already-known solution (lowerLimit), one can
 * immediately prune the decision tree --- discard this partial solution.
 *
 * Branch-and-Bound solution, using a "best-fit-first" rule for expanding
 * states.  As a maximization problem, it has been verified by experiment
 * that the priority queue should be driven by the present fixed value of
 * the permutation, not the upper-limit value discussed above.
 *
 * Language:  Java, version 5 or later (Scanner, printf)
 *
 * Author:  Timothy Rolfe
 */

#include <time.h>       // clock() etc.
#include <stdio.h>
#include <string.h>     // memcpy()
#include <stdlib.h>
#include "State.h"
#include "MinHeap.h"

#ifdef _WIN32           // NOT the Unix environment
// Use a macro to get double seconds from clock()
#define cpuClock()  ((double)clock()/CLOCKS_PER_SEC)
#define wallClock() ((double)clock()/CLOCKS_PER_SEC)
#elif defined(__TURBOC__)
#define cpuClock()  ((double)clock()/CLOCKS_PER_SEC)
#define wallClock() ((double)clock()/CLOCKS_PER_SEC)
#else
#include "cpuTimes.c"
#endif

// Diagonal/antidiagonal permutation sets solution and lowerLimit
#define MAX_SIZE 35
int solution[MAX_SIZE], lowerLimit;
int size, benefit[MAX_SIZE][MAX_SIZE];
enum { FALSE, TRUE };
int TRACE = FALSE;

void initialize(FILE *fptr)
{
   int row, col, k,
       sum1 = 0, sum2 = 0;

// benefit and solution are global fixed-size arrays
   fscanf(fptr, "%d", &size);

   for (row = 0 , k = size; row < size; row++)
   {  solution[row] = row;
      for (col = 0; col < size; col++)
         fscanf (fptr, "%d", &benefit[row][col]);
      sum1 += benefit[row][row];
      sum2 += benefit[row][--k];     // Note predecrement of k
   }
   if (sum1 >= sum2)
      lowerLimit = sum1;
   else
   {  lowerLimit = sum2;
      col = size;
      for (row = 0; row < size; row++)
         solution[row] = --col;  // Antidiagonal
   }
}

static void displayCost()
{
   int row, col;

   printf("Size %d\n", size);
   for (row = 0; row < size; row++)
   {  for (col = 0; col < size; col++)
         printf("%4d", benefit[row][col]);
      putchar('\n');
   }
   printf("Lower bound: %d\nVector: ", lowerLimit);
   for (row = 0; row < size; row++)
      printf("%3d", solution[row]);
   putchar('\n');
}

void swap (int *x, int p, int q)
{  int t = x[p]; x[p] = x[q]; x[q] = t;  }

int colMaxSum(int start, int *work)
{  int sum = 0;
   for (int k = start; k < size; k++)
   {
      int row, col = work[k],
          columnMaximum = benefit[start][col];
      for (row = start+1; row < size; row++)
         if (columnMaximum < benefit[row][col])
             columnMaximum = benefit[row][col];
      sum += columnMaximum;
   }
   return sum;
}

// Working from the global variables, obtain a maximal-benefit solution.
// The initialize method has already computing the initial upper-limit
// as the diagonal or antidiagonal permutation
void explore()
{
   MinHeap pq(size);
   int    *work = (int*) calloc (size, sizeof *work);

   pq.add( new State(solution, size, 0, 0, colMaxSum(0, solution)));
   while (!pq.isEmpty())
   {
      State* eNode = pq.remove();
      int    score = eNode->_score,
             myMax = eNode->_myMax,
             idx   = eNode->_idx;
      int    j;

      memcpy (work, eNode->_vector, size * sizeof *work);
      if (TRACE)
      {  char *pict = eNode->toString();
         printf("Expanding %s\n", pict);
         free(pict);
      }
      if (myMax <= lowerLimit)
      {  if (TRACE) printf("Lower limit %d, so discard\n",
             lowerLimit);
         continue;      // I.e., discard
      }
      // Permute available scores through [idx]
      for ( j = idx; j < size; j++ )
      {  int workScore = score;
         State* addendum;

         swap (work, idx, j);
         workScore += benefit[idx][work[idx]];
         // Compute the column maxima to the right of idx
         myMax = workScore + colMaxSum(idx+1, work);
         // Selecting [size-2] also fixes [size-1] ---
         // a permutation has been completed.
         if ( idx == size-2 )
         {
            if (myMax > lowerLimit)
            {
            // System.arraycopy(work, 0, solution, 0, size);
               memcpy(solution, work, size*sizeof *work);
               if (TRACE) printf("Accept, score > %d\n",
                                            lowerLimit);
               lowerLimit = myMax;
            }
            else if (TRACE) printf("Reject, score <= %d\n",
                                              lowerLimit);
         }
         else
         {
            addendum = new State(work, size, idx+1, workScore, myMax);
            if (myMax > lowerLimit)
            {
                pq.add(addendum);
                if (TRACE)
                   printf("\tMove forward with %s\n", addendum->toString());
             }
            else if (TRACE)
            {  printf("\tLower bound %d, so reject %s\n",
                  lowerLimit, addendum->toString());
               delete addendum;
            }
         }
      }
      delete eNode;     // Avoid a memory leak
   }
// if (TRACE)
      printf ("Final size of the MinHeap:  %d\n", pq.size()); fflush(stdout);
   // Work vector returned to the heap
   free (work);
}

static void displaySolution()
{
   int idx;

   printf("Branch-and-Bound sequential: %d\n", lowerLimit);
   for (idx = 0; idx < size; idx++)
      printf ("%3d", solution[idx]);
   putchar('\n');
}

int main(int argc, char* argv[])
{
   const char *filename = argc < 2 ? "Asg5.in" : argv[1];
   FILE   *fptr = fopen(filename, "r");
   double start, elapsed;

   printf("Reading from %s\n", filename);
   if (fptr == NULL)
   {  puts("Failed to open file."); return 127;  }

   initialize(fptr);
   if (TRACE)
      displayCost();
   start = cpuClock();
   explore();
   elapsed = cpuClock() - start;
   displaySolution();
   printf("Time:  %3.3f seconds\n", elapsed);
   return 0;
}