/*
 * ex1_knapsack.c -- 0/1 knapsack
 *
 *  max   sum  v_i x_i
 *  s.t.  sum  w_i x_i  <=  W
 *        x_i in {0, 1}
 *
 * What the LP relaxation buys you:
 *
 *   The LP solution to a 0/1 knapsack is Dantzig's classical fractional
 *   knapsack: sort items by v_i/w_i, fill greedily, and the last item
 *   that doesn't fit goes in fractionally. This is provably the tightest
 *   linear upper bound on the integer optimum -- there is at most one
 *   fractional variable, every other variable is already 0 or 1, and the
 *   gap to the integer optimum is bounded by the value of the most
 *   valuable single item.
 *
 *   Branch-and-bound therefore has very little work to do: each node's
 *   LP either rounds the one fractional variable down (drop the item)
 *   or up (force-include it), and one of the two children has an LP
 *   bound at most v_i/w_i below the parent. The tree is shallow and
 *   most subtrees prune immediately against the incumbent.
 */
#include "../src/milp.h"
#include <stdio.h>
#include <stdlib.h>

int main(void) {
    /* 8 items, capacity W = 12.
     * Designed so the LP rounds at exactly one fractional variable AND the
     * integer optimum sits strictly below the LP bound -- so branch-and-bound
     * actually has to explore both branches of the fractional variable. */
    enum { N = 8 };
    double values [N] = { 10, 8, 7, 6, 5, 4, 3, 2 };
    double weights[N] = {  5, 4, 4, 4, 3, 3, 2, 1 };
    double W = 12;

    /* Decision: x_i in {0,1} for each item.
     * Single capacity row. */
    MilpModel *m = milp_new(N, 1);
    milp_set_objective(m, MILP_MAX, values);
    milp_set_row(m, 0, weights, MILP_LE, W);
    for (int j = 0; j < N; j++) milp_set_var_type(m, j, MILP_BINARY);

    /* --- LP relaxation --- */
    double x_lp[N], obj_lp;
    milp_solve_lp(m, x_lp, &obj_lp);
    printf("=== LP relaxation ===\n");
    printf("  objective = %.4f\n", obj_lp);
    printf("  x = [");
    for (int j = 0; j < N; j++) printf(" %.3f", x_lp[j]);
    printf(" ]\n");
    int n_frac = 0;
    for (int j = 0; j < N; j++) {
        double f = x_lp[j] - (int)x_lp[j];
        if (f > 1e-6 && f < 1 - 1e-6) n_frac++;
    }
    printf("  fractional variables: %d  (Dantzig's bound: at most 1)\n\n", n_frac);

    /* --- MIP --- */
    double x[N], obj;
    int s = milp_solve(m, x, &obj);
    printf("=== Integer solution ===\n");
    printf("  status     = %d\n", s);
    printf("  objective  = %.4f\n", obj);
    printf("  items in   = ");
    double w_used = 0;
    for (int j = 0; j < N; j++) {
        if (x[j] > 0.5) {
            printf("%d ", j);
            w_used += weights[j];
        }
    }
    printf("(weight %.0f / %.0f)\n\n", w_used, W);

    printf("LP bound       %.4f\n", milp_lp_relaxation(m));
    printf("MIP optimum    %.4f\n", obj);
    printf("Integrality gap %.4f  (%.2f%%)\n",
           milp_lp_relaxation(m) - obj,
           100.0 * (milp_lp_relaxation(m) - obj) /
           (obj > 1e-9 ? obj : 1));
    printf("B&B nodes      %d\n", milp_nodes_explored(m));
    printf("LPs solved     %d\n", milp_lps_solved(m));

    milp_free(m);
    return 0;
}
