/*
 * ex3_setcover.c -- minimum set cover
 *
 *  min   sum  x_s
 *  s.t.  sum_{s : e in S_s}  x_s  >=  1   for each element e
 *        x_s in {0, 1}
 *
 * What the LP relaxation buys you:
 *
 *   Set cover is NP-hard, and unlike assignment it has a real integrality
 *   gap. The LP relaxation produces a *fractional cover* in which a set
 *   may be "half chosen", which is pleasant to think about but useless to
 *   ship. Yet the LP value is still a valid lower bound on the integer
 *   optimum, and a famous result (Lovasz, Chvatal) tells us that rounding
 *   each x_s up to 1 wherever its LP value exceeds 1/f -- with f the
 *   maximum element frequency -- gives a cover within an O(log n) factor
 *   of optimum.
 *
 *   For us the LP serves two distinct purposes: it is the bound that
 *   prunes the branch-and-bound tree, and it is a starting heuristic
 *   that often produces a feasible incumbent before any branching.
 *   This example shows both: the LP value, the rounded heuristic, and
 *   the optimum reached by branch-and-bound.
 */
#include "../src/milp.h"
#include <stdio.h>
#include <stdlib.h>

int main(void) {
    /* Universe of 8 elements covered by 6 sets. */
    enum { N_ELEMS = 8, N_SETS = 6 };
    int sets[N_SETS][8] = {
        /* S0 */ { 0, 1, 2, 3, -1 },
        /* S1 */ { 0, 4, 5,    -1 },
        /* S2 */ { 1, 4, 6,    -1 },
        /* S3 */ { 2, 5, 7,    -1 },
        /* S4 */ { 3, 6, 7,    -1 },
        /* S5 */ { 0, 1, 4, 5, 6, -1 },
    };

    /* One row per element: sum of x_s over sets containing e is >= 1. */
    MilpModel *m = milp_new(N_SETS, N_ELEMS);

    double c[N_SETS];
    for (int s = 0; s < N_SETS; s++) c[s] = 1.0;   /* unit costs */
    milp_set_objective(m, MILP_MIN, c);

    for (int e = 0; e < N_ELEMS; e++) {
        double row[N_SETS];
        for (int s = 0; s < N_SETS; s++) row[s] = 0.0;
        for (int s = 0; s < N_SETS; s++)
            for (int k = 0; sets[s][k] >= 0; k++)
                if (sets[s][k] == e) row[s] = 1.0;
        milp_set_row(m, e, row, MILP_GE, 1.0);
    }

    for (int s = 0; s < N_SETS; s++) milp_set_var_type(m, s, MILP_BINARY);

    /* --- LP relaxation --- */
    double x_lp[N_SETS], obj_lp;
    milp_solve_lp(m, x_lp, &obj_lp);
    printf("=== LP relaxation ===\n");
    printf("  objective = %.4f  (fractional lower bound)\n", obj_lp);
    printf("  x =");
    for (int s = 0; s < N_SETS; s++) printf(" %.3f", x_lp[s]);
    printf("\n\n");

    /* --- LP rounding heuristic: include S_s if x_s >= 1/f.
     * f = max element frequency = how many sets cover the most-covered element. */
    int max_freq = 0;
    for (int e = 0; e < N_ELEMS; e++) {
        int freq = 0;
        for (int s = 0; s < N_SETS; s++)
            for (int k = 0; sets[s][k] >= 0; k++)
                if (sets[s][k] == e) { freq++; break; }
        if (freq > max_freq) max_freq = freq;
    }
    double thresh = 1.0 / max_freq;
    int rounded = 0;
    printf("=== LP rounding heuristic (threshold 1/%d = %.3f) ===\n",
           max_freq, thresh);
    printf("  picked: ");
    for (int s = 0; s < N_SETS; s++) {
        if (x_lp[s] >= thresh - 1e-9) {
            printf("S%d ", s);
            rounded++;
        }
    }
    printf("\n  size   = %d  (guaranteed within factor %d of optimum)\n\n",
           rounded, max_freq);

    /* --- Exact MIP --- */
    double x[N_SETS], obj;
    milp_solve(m, x, &obj);
    printf("=== Integer optimum (branch-and-bound) ===\n");
    printf("  objective = %.4f\n", obj);
    printf("  picked   :");
    for (int s = 0; s < N_SETS; s++) if (x[s] > 0.5) printf(" S%d", s);
    printf("\n\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\n", milp_nodes_explored(m));
    printf("LPs solved     %d\n", milp_lps_solved(m));

    milp_free(m);
    return 0;
}
