/* sanity.c -- quick smoke-test of milp.c */
#include "src/milp.h"
#include <stdio.h>
#include <math.h>
#include <stdlib.h>

static int approx(double a, double b) { return fabs(a - b) < 1e-5; }

static int test_lp(void) {
    /* max x + y  s.t.  x + 2y <= 4,  2x + y <= 4,  x,y >= 0
     * Optimal: x=y=4/3, obj=8/3. */
    MilpModel *m = milp_new(2, 2);
    double c[] = {1, 1};
    milp_set_objective(m, MILP_MAX, c);
    double r1[] = {1, 2}; milp_set_row(m, 0, r1, MILP_LE, 4);
    double r2[] = {2, 1}; milp_set_row(m, 1, r2, MILP_LE, 4);

    double x[2], obj;
    int s = milp_solve_lp(m, x, &obj);
    printf("LP test 1: status=%d, obj=%.6f, x=[%.6f, %.6f]\n", s, obj, x[0], x[1]);
    int ok = (s == MILP_OPTIMAL && approx(obj, 8.0 / 3.0) &&
              approx(x[0], 4.0 / 3.0) && approx(x[1], 4.0 / 3.0));
    milp_free(m);
    return ok;
}

static int test_min(void) {
    /* min 3x + 2y  s.t.  x + y >= 4, x,y >= 0
     * Optimal: x=0, y=4, obj=8. */
    MilpModel *m = milp_new(2, 1);
    double c[] = {3, 2};
    milp_set_objective(m, MILP_MIN, c);
    double r1[] = {1, 1}; milp_set_row(m, 0, r1, MILP_GE, 4);

    double x[2], obj;
    int s = milp_solve_lp(m, x, &obj);
    printf("LP test 2: status=%d, obj=%.6f, x=[%.6f, %.6f]\n", s, obj, x[0], x[1]);
    int ok = (s == MILP_OPTIMAL && approx(obj, 8.0));
    milp_free(m);
    return ok;
}

static int test_eq(void) {
    /* min x + y  s.t.  x + y == 5,  x - y >= 1
     * Optimal: x=3, y=2, obj=5. */
    MilpModel *m = milp_new(2, 2);
    double c[] = {1, 1};
    milp_set_objective(m, MILP_MIN, c);
    double r1[] = {1, 1};  milp_set_row(m, 0, r1, MILP_EQ, 5);
    double r2[] = {1, -1}; milp_set_row(m, 1, r2, MILP_GE, 1);

    double x[2], obj;
    int s = milp_solve_lp(m, x, &obj);
    printf("LP test 3: status=%d, obj=%.6f, x=[%.6f, %.6f]\n", s, obj, x[0], x[1]);
    int ok = (s == MILP_OPTIMAL && approx(obj, 5.0) && approx(x[0] + x[1], 5.0));
    milp_free(m);
    return ok;
}

static int test_mip(void) {
    /* max 5x + 4y  s.t.  6x + 4y <= 24, x + 2y <= 6, x,y in {0,1,2,3,4}
     * LP optimum: x=3, y=1.5 -> obj 21
     * MIP optimum: x=3, y=1 (or x=2, y=2) -> obj 19, depending on tie-breaking,
     * actually x=2, y=2 gives 18; x=3, y=1 gives 19; x=4, y=0 gives 20.
     * Wait: 6*4 + 4*0 = 24 OK, 1*4 + 2*0 = 4 OK. So x=4, y=0 -> obj 20.
     * Actually re-checking: x=3, y=1: 6*3+4=22<=24, 3+2=5<=6, obj=19.
     * x=4, y=0: 24<=24, 4<=6, obj=20. So integer optimum is 20.
     * (Or x=3, y=1.5 LP gives 21, then branch.) */
    MilpModel *m = milp_new(2, 2);
    double c[] = {5, 4};
    milp_set_objective(m, MILP_MAX, c);
    double r1[] = {6, 4}; milp_set_row(m, 0, r1, MILP_LE, 24);
    double r2[] = {1, 2}; milp_set_row(m, 1, r2, MILP_LE, 6);
    milp_set_var_type(m, 0, MILP_INTEGER);
    milp_set_var_type(m, 1, MILP_INTEGER);
    milp_set_var_bounds(m, 0, 0, 4);
    milp_set_var_bounds(m, 1, 0, 4);

    double x[2], obj;
    int s = milp_solve(m, x, &obj);
    printf("MIP test:  status=%d, obj=%.6f, x=[%.6f, %.6f], LP=%.6f, nodes=%d\n",
           s, obj, x[0], x[1], milp_lp_relaxation(m), milp_nodes_explored(m));
    int ok = (s == MILP_OPTIMAL && approx(obj, 20.0));
    milp_free(m);
    return ok;
}

static int test_infeasible(void) {
    /* x + y <= 1, x + y >= 3  -> infeasible */
    MilpModel *m = milp_new(2, 2);
    double c[] = {1, 1};
    milp_set_objective(m, MILP_MIN, c);
    double r1[] = {1, 1}; milp_set_row(m, 0, r1, MILP_LE, 1);
    double r2[] = {1, 1}; milp_set_row(m, 1, r2, MILP_GE, 3);

    double x[2], obj;
    int s = milp_solve_lp(m, x, &obj);
    printf("Infeasible test: status=%d (expect %d)\n", s, MILP_INFEASIBLE);
    milp_free(m);
    return s == MILP_INFEASIBLE;
}

int main(void) {
    int p = 0, n = 0;
    n++; if (test_lp())         p++; else printf("  FAILED\n");
    n++; if (test_min())        p++; else printf("  FAILED\n");
    n++; if (test_eq())         p++; else printf("  FAILED\n");
    n++; if (test_mip())        p++; else printf("  FAILED\n");
    n++; if (test_infeasible()) p++; else printf("  FAILED\n");
    printf("\n%d / %d passed\n", p, n);
    return p == n ? 0 : 1;
}
