/*
 * ex4_wta.c -- budgeted weapon-target assignment (linear WTA)
 *
 *  max   sum_{i,j}  v_j p_ij x_ij                  (expected damage)
 *  s.t.  sum_j  x_ij        <=  1     each weapon engages at most one target
 *        sum_i  x_ij        <=  1     each target engaged by at most one weapon
 *        sum_{i,j}  c_ij x_ij  <=  B  total cost (rounds, fuel, anything) <= B
 *        x_ij in {0, 1}
 *
 * What the LP relaxation buys you:
 *
 *   Without the budget row, the constraint matrix is the bipartite
 *   incidence matrix -- totally unimodular, just like ex2. Add a single
 *   "knapsack" budget row and TUM is destroyed; the LP can return
 *   fractional engagements that "share" the budget across targets in a
 *   way no real shooter can. The LP optimum sits strictly above the MIP
 *   optimum, but the gap is small because only the budget row is loose.
 *   Branch-and-bound prunes most of the 2^{n*m} space using this bound.
 *
 *   This is the kernel that sits inside a static-WTA solver: in a real
 *   TEWA system one wraps it in a time loop, recomputes p_ij from
 *   geometry and seeker models, and re-solves at each engagement update.
 */
#include "../src/milp.h"
#include <stdio.h>
#include <stdlib.h>

int main(void) {
    enum { NW = 4, NT = 3 };           /* weapons, targets */
    double v[NT] = { 10.0, 8.0, 6.0 }; /* target values */

    /* Kill probability p[i][j] for weapon i firing at target j. */
    double p[NW][NT] = {
        { 0.9, 0.7, 0.5 },
        { 0.8, 0.9, 0.4 },
        { 0.6, 0.5, 0.9 },
        { 0.7, 0.6, 0.5 },
    };

    /* Cost (rounds, fuel, time -- whatever the doctrine says). */
    double cost[NW][NT] = {
        { 3, 3, 2 },
        { 3, 4, 2 },
        { 2, 2, 3 },
        { 2, 2, 1 },
    };
    double B = 7.0;

    /* Variables: x_{i,j} indexed as i*NT + j. */
    int n_vars = NW * NT;
    int n_rows = NW + NT + 1;          /* weapon caps + target caps + budget */
    MilpModel *m = milp_new(n_vars, n_rows);

    /* Objective: maximise expected damage. */
    double *c = malloc(sizeof(double) * n_vars);
    for (int i = 0; i < NW; i++)
        for (int j = 0; j < NT; j++)
            c[i * NT + j] = v[j] * p[i][j];
    milp_set_objective(m, MILP_MAX, c);

    /* Each weapon used at most once. */
    double *row = calloc(n_vars, sizeof(double));
    int r = 0;
    for (int i = 0; i < NW; i++) {
        for (int k = 0; k < n_vars; k++) row[k] = 0.0;
        for (int j = 0; j < NT; j++) row[i * NT + j] = 1.0;
        milp_set_row(m, r++, row, MILP_LE, 1.0);
    }
    /* Each target hit at most once. */
    for (int j = 0; j < NT; j++) {
        for (int k = 0; k < n_vars; k++) row[k] = 0.0;
        for (int i = 0; i < NW; i++) row[i * NT + j] = 1.0;
        milp_set_row(m, r++, row, MILP_LE, 1.0);
    }
    /* Budget. */
    for (int k = 0; k < n_vars; k++) row[k] = 0.0;
    for (int i = 0; i < NW; i++)
        for (int j = 0; j < NT; j++)
            row[i * NT + j] = cost[i][j];
    milp_set_row(m, r++, row, MILP_LE, B);
    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);

    printf("=== LP relaxation ===\n");
    printf("  expected damage (UB) = %.4f\n", obj_lp);
    printf("  fractional engagements:\n");
    int n_frac = 0;
    for (int i = 0; i < NW; i++) {
        for (int j = 0; j < NT; j++) {
            double x = x_lp[i * NT + j];
            if (x > 1e-6 && x < 1 - 1e-6) {
                printf("    W%d -> T%d : %.3f\n", i, j, x);
                n_frac++;
            }
        }
    }
    if (n_frac == 0) printf("    (none)\n");
    printf("\n");

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

    printf("=== Integer engagement plan ===\n");
    printf("  expected damage = %.4f\n", obj);
    printf("  fire orders:\n");
    double cost_used = 0;
    for (int i = 0; i < NW; i++) {
        for (int j = 0; j < NT; j++) {
            if (x[i * NT + j] > 0.5) {
                printf("    W%d -> T%d   p=%.2f  v=%.0f  cost=%.0f\n",
                       i, j, p[i][j], v[j], cost[i][j]);
                cost_used += cost[i][j];
            }
        }
    }
    printf("  cost spent  = %.0f / %.0f\n\n", cost_used, B);

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

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