/*
 * ex2_assignment.c -- linear assignment (n agents, n tasks)
 *
 *  min   sum  c_ij x_ij
 *  s.t.  sum_j  x_ij  =  1     for each agent i
 *        sum_i  x_ij  =  1     for each task  j
 *        x_ij in {0, 1}
 *
 * What the LP relaxation buys you:
 *
 *   The constraint matrix of an n x n assignment problem is *totally
 *   unimodular*: every square submatrix has determinant 0, +1, or -1.
 *   A theorem of Hoffman & Kruskal then tells us that for any integer
 *   right-hand side (here, all 1s), every basic feasible LP solution is
 *   already integer. The LP relaxation is exact -- no branching is ever
 *   needed.
 *
 *   This is the strongest possible relaxation: gap = 0. Branch-and-bound
 *   solves the problem at the root with one LP and zero branches. In
 *   production you would use the Hungarian or auction algorithm, which
 *   exploit the structure directly; the LP simply confirms that the
 *   problem is "easy" in a deep sense.
 */
#include "../src/milp.h"
#include <stdio.h>
#include <stdlib.h>

int main(void) {
    /* 4 agents, 4 tasks. cost[i][j] = c_ij. */
    enum { N = 4 };
    double cost[N][N] = {
        {  9,  2,  7,  8 },
        {  6,  4,  3,  7 },
        {  5,  8,  1,  8 },
        {  7,  6,  9,  4 },
    };

    /* Variables: x_{i,j} indexed as i*N + j. */
    int n_vars = N * N;
    int n_rows = 2 * N;          /* N row sums + N column sums */
    MilpModel *m = milp_new(n_vars, n_rows);

    double *c = malloc(sizeof(double) * n_vars);
    for (int i = 0; i < N; i++)
        for (int j = 0; j < N; j++)
            c[i * N + j] = cost[i][j];
    milp_set_objective(m, MILP_MIN, c);

    /* Each agent does exactly one task. */
    double *row = calloc(n_vars, sizeof(double));
    for (int i = 0; i < N; i++) {
        for (int k = 0; k < n_vars; k++) row[k] = 0.0;
        for (int j = 0; j < N; j++) row[i * N + j] = 1.0;
        milp_set_row(m, i, row, MILP_EQ, 1.0);
    }
    /* Each task done by exactly one agent. */
    for (int j = 0; j < N; j++) {
        for (int k = 0; k < n_vars; k++) row[k] = 0.0;
        for (int i = 0; i < N; i++) row[i * N + j] = 1.0;
        milp_set_row(m, N + j, row, MILP_EQ, 1.0);
    }
    free(row);
    free(c);

    for (int j = 0; j < n_vars; j++) milp_set_var_type(m, j, MILP_BINARY);

    /* --- LP relaxation --- */
    double *x_lp = malloc(sizeof(double) * n_vars);
    double obj_lp;
    milp_solve_lp(m, x_lp, &obj_lp);

    int n_frac = 0;
    for (int j = 0; j < n_vars; j++) {
        double f = x_lp[j] - (int)x_lp[j];
        if (f > 1e-6 && f < 1 - 1e-6) n_frac++;
    }
    printf("=== LP relaxation ===\n");
    printf("  objective           = %.4f\n", obj_lp);
    printf("  fractional variables = %d   (TUM predicts 0)\n\n", n_frac);

    /* --- MIP --- */
    double *x = malloc(sizeof(double) * n_vars);
    double obj;
    milp_solve(m, x, &obj);

    printf("=== Integer solution ===\n");
    printf("  objective = %.4f\n", obj);
    printf("  assignment:\n");
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            if (x[i * N + j] > 0.5)
                printf("    agent %d -> task %d  (cost %.0f)\n",
                       i, j, cost[i][j]);
        }
    }
    printf("\n");
    printf("LP bound       %.4f\n", milp_lp_relaxation(m));
    printf("MIP optimum    %.4f\n", obj);
    printf("Integrality gap %.4f\n", obj - milp_lp_relaxation(m));
    printf("B&B nodes      %d   (1 = LP was already integer)\n",
           milp_nodes_explored(m));
    printf("LPs solved     %d\n", milp_lps_solved(m));

    free(x_lp);
    free(x);
    milp_free(m);
    return 0;
}
