pixie.c

/*
Description:
PIXIE is a compact mixed-integer linear programming solver implemented in one ISO C17 file.
It includes an LP/MPS reader, a two-phase primal simplex core, and a deterministic branch-and-bound layer for integer and binary variables.
This port preserves the original solver behavior while removing non-standard dependencies and keeping the implementation portable and self-contained.

Date: 2026-02-06
Author: Riveros / ASPIE / SATX
Copyright: © Oscar Riveros. All rights reserved.

gcc -O3 -std=c17 -Wall -Wextra -pedantic pixie.c -o pixie
*/

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdbool.h>
#include <string.h>
#include <math.h>
#include <ctype.h>
#include <errno.h>
#include <limits.h>
#include <time.h>

#define PIXIE_FEAS_TOL 1e-8
#define PIXIE_INT_TOL 1e-7
#define PIXIE_PIV_TOL 1e-10
#define PIXIE_EPS 1e-12

typedef enum {
    PIXIE_CMP_LE = -1,
    PIXIE_CMP_EQ = 0,
    PIXIE_CMP_GE = 1
} PixieCompare;

typedef enum {
    PIXIE_VAR_CONT = 0,
    PIXIE_VAR_INT = 1,
    PIXIE_VAR_BIN = 2
} PixieVarType;

typedef enum {
    PIXIE_STATUS_OPTIMAL = 0,
    PIXIE_STATUS_INFEASIBLE = 1,
    PIXIE_STATUS_UNBOUNDED = 2,
    PIXIE_STATUS_UNKNOWN = 3,
    PIXIE_STATUS_ERROR = 4
} PixieStatus;

typedef enum {
    PIXIE_FORMAT_AUTO = 0,
    PIXIE_FORMAT_LP = 1,
    PIXIE_FORMAT_MPS = 2
} PixieFormat;

typedef struct {
    int *idx;
    double *val;
    int len;
    int cap;
} PixieSparse;

typedef struct {
    PixieSparse a;
    PixieCompare cmp;
    double rhs;
} PixieConstraint;

typedef struct {
    char *name;
    double obj;
    double lb;
    double ub;
    PixieVarType type;
} PixieVar;

typedef struct {
    PixieVar *vars;
    int n_vars;
    int cap_vars;
    PixieConstraint *cons;
    int n_cons;
    int cap_cons;
    int obj_sense; /* +1 minimize, -1 maximize */
} PixieModel;

typedef struct {
    PixieFormat format;
    const char *file;
    double time_limit_sec;   /* <= 0 means unlimited */
    long long node_limit;    /* <= 0 means unlimited */
    double gap_limit;        /* < 0 means disabled */
    uint64_t seed;
    bool seed_set;
    int verbose;             /* 0..3 */
    bool pure_lp;
    bool selftest;
} PixieOptions;

typedef struct {
    PixieStatus status;
    double obj_min;
    double obj_out;
    double *x;
    int n;
    bool has_primal;
    long long nodes_processed;
    bool stopped_time;
    bool stopped_nodes;
    bool stopped_gap;
} PixieSolution;

typedef struct {
    double **table;
    double *raw;
    int *v_idx;
    bool *a_flg;
    int v_cnt;
    int e_cnt;
    int s_cnt;
    int a_cnt;
} PixieSimplex;

typedef enum {
    PIXIE_SIMPLEX_OPTIMAL = 0,
    PIXIE_SIMPLEX_UNBOUNDED = 1,
    PIXIE_SIMPLEX_NUMERIC = 2
} PixieSimplexStatus;

typedef struct {
    PixieConstraint *data;
    int len;
    int cap;
} PixieConVec;

typedef struct {
    double *lb;
    double *ub;
    double bound;
    int depth;
} PixieNode;

typedef struct {
    PixieNode *data;
    int len;
    int cap;
} PixieNodeStack;

typedef struct {
    char **data;
    int len;
    int cap;
} PixieTokVec;

typedef enum {
    PIXIE_MAP_SHIFT_LB = 0,
    PIXIE_MAP_SHIFT_UB = 1,
    PIXIE_MAP_FREE = 2
} PixieVarMapKind;

typedef struct {
    PixieVarMapKind kind;
    int y0;
    int y1;
    double constant;
} PixieVarMap;

static void pixie_set_error(char *err, size_t errsz, const char *msg) {
    if (err != NULL && errsz > 0) {
        size_t n = strlen(msg);
        if (n >= errsz) {
            n = errsz - 1;
        }
        memcpy(err, msg, n);
        err[n] = '\0';
    }
}

static void *pixie_xmalloc(size_t n) {
    void *p = malloc(n == 0 ? 1 : n);
    if (p == NULL) {
        fprintf(stderr, "fatal: out of memory\n");
        exit(EXIT_FAILURE);
    }
    return p;
}

static void *pixie_xcalloc(size_t count, size_t size) {
    void *p = calloc(count == 0 ? 1 : count, size == 0 ? 1 : size);
    if (p == NULL) {
        fprintf(stderr, "fatal: out of memory\n");
        exit(EXIT_FAILURE);
    }
    return p;
}

static void *pixie_xrealloc(void *ptr, size_t n) {
    void *p = realloc(ptr, n == 0 ? 1 : n);
    if (p == NULL) {
        fprintf(stderr, "fatal: out of memory\n");
        exit(EXIT_FAILURE);
    }
    return p;
}

static char *pixie_strdup_c(const char *s) {
    size_t n = strlen(s);
    char *d = (char *)pixie_xmalloc(n + 1);
    memcpy(d, s, n + 1);
    return d;
}

static int pixie_stricmp_c(const char *a, const char *b) {
    unsigned char ca;
    unsigned char cb;
    while (*a != '\0' && *b != '\0') {
        ca = (unsigned char)tolower((unsigned char)*a);
        cb = (unsigned char)tolower((unsigned char)*b);
        if (ca != cb) {
            return (int)ca - (int)cb;
        }
        ++a;
        ++b;
    }
    ca = (unsigned char)tolower((unsigned char)*a);
    cb = (unsigned char)tolower((unsigned char)*b);
    return (int)ca - (int)cb;
}

static bool pixie_starts_with_ci(const char *s, const char *prefix) {
    while (*prefix != '\0') {
        if (*s == '\0') {
            return false;
        }
        if (tolower((unsigned char)*s) != tolower((unsigned char)*prefix)) {
            return false;
        }
        ++s;
        ++prefix;
    }
    return true;
}

static bool pixie_ends_with_ci(const char *s, const char *suffix) {
    size_t n = strlen(s);
    size_t m = strlen(suffix);
    if (m > n) {
        return false;
    }
    return pixie_stricmp_c(s + (n - m), suffix) == 0;
}

static char *pixie_ltrim(char *s) {
    while (*s != '\0' && isspace((unsigned char)*s)) {
        ++s;
    }
    return s;
}

static void pixie_rtrim(char *s) {
    size_t n = strlen(s);
    while (n > 0) {
        if (!isspace((unsigned char)s[n - 1])) {
            break;
        }
        s[n - 1] = '\0';
        --n;
    }
}

static char *pixie_trim(char *s) {
    char *t = pixie_ltrim(s);
    pixie_rtrim(t);
    return t;
}

static bool pixie_is_finite(double x) {
    return isfinite(x) != 0;
}

static double pixie_max(double a, double b) {
    return (a > b) ? a : b;
}

static double pixie_min(double a, double b) {
    return (a < b) ? a : b;
}

static bool pixie_parse_double_token(const char *s, double *out) {
    char *end = NULL;
    double v;
    errno = 0;
    v = strtod(s, &end);
    if (end == s) {
        return false;
    }
    while (*end != '\0' && isspace((unsigned char)*end)) {
        ++end;
    }
    if (*end != '\0') {
        return false;
    }
    if (errno == ERANGE) {
        return false;
    }
    *out = v;
    return true;
}

static bool pixie_parse_ll_token(const char *s, long long *out) {
    char *end = NULL;
    long long v;
    errno = 0;
    v = strtoll(s, &end, 10);
    if (end == s) {
        return false;
    }
    while (*end != '\0' && isspace((unsigned char)*end)) {
        ++end;
    }
    if (*end != '\0') {
        return false;
    }
    if (errno == ERANGE) {
        return false;
    }
    *out = v;
    return true;
}

static bool pixie_parse_u64_token(const char *s, uint64_t *out) {
    char *end = NULL;
    unsigned long long v;
    errno = 0;
    v = strtoull(s, &end, 10);
    if (end == s) {
        return false;
    }
    while (*end != '\0' && isspace((unsigned char)*end)) {
        ++end;
    }
    if (*end != '\0') {
        return false;
    }
    if (errno == ERANGE) {
        return false;
    }
    *out = (uint64_t)v;
    return true;
}

static uint64_t pixie_rng_next(uint64_t *state) {
    uint64_t x = *state;
    if (x == 0) {
        x = UINT64_C(0x9e3779b97f4a7c15);
    }
    x ^= x >> 12;
    x ^= x << 25;
    x ^= x >> 27;
    *state = x;
    return x * UINT64_C(2685821657736338717);
}

static void pixie_sparse_init(PixieSparse *s) {
    s->idx = NULL;
    s->val = NULL;
    s->len = 0;
    s->cap = 0;
}

static void pixie_sparse_free(PixieSparse *s) {
    free(s->idx);
    free(s->val);
    s->idx = NULL;
    s->val = NULL;
    s->len = 0;
    s->cap = 0;
}

static void pixie_sparse_reserve(PixieSparse *s, int need) {
    int nc;
    if (need <= s->cap) {
        return;
    }
    nc = (s->cap <= 0) ? 8 : s->cap;
    while (nc < need) {
        if (nc > INT_MAX / 2) {
            nc = need;
            break;
        }
        nc *= 2;
    }
    s->idx = (int *)pixie_xrealloc(s->idx, (size_t)nc * sizeof(int));
    s->val = (double *)pixie_xrealloc(s->val, (size_t)nc * sizeof(double));
    s->cap = nc;
}

static void pixie_sparse_add(PixieSparse *s, int idx, double val) {
    int i;
    if (fabs(val) <= PIXIE_EPS) {
        return;
    }
    for (i = 0; i < s->len; ++i) {
        if (s->idx[i] == idx) {
            s->val[i] += val;
            return;
        }
    }
    pixie_sparse_reserve(s, s->len + 1);
    s->idx[s->len] = idx;
    s->val[s->len] = val;
    ++s->len;
}

static void pixie_sparse_prune(PixieSparse *s, double tol) {
    int i;
    int w = 0;
    for (i = 0; i < s->len; ++i) {
        if (fabs(s->val[i]) > tol) {
            s->idx[w] = s->idx[i];
            s->val[w] = s->val[i];
            ++w;
        }
    }
    s->len = w;
}

static void pixie_sparse_move(PixieSparse *dst, PixieSparse *src) {
    dst->idx = src->idx;
    dst->val = src->val;
    dst->len = src->len;
    dst->cap = src->cap;
    src->idx = NULL;
    src->val = NULL;
    src->len = 0;
    src->cap = 0;
}

static void pixie_sparse_clone(PixieSparse *dst, const PixieSparse *src) {
    int i;
    pixie_sparse_init(dst);
    if (src->len <= 0) {
        return;
    }
    pixie_sparse_reserve(dst, src->len);
    for (i = 0; i < src->len; ++i) {
        dst->idx[i] = src->idx[i];
        dst->val[i] = src->val[i];
    }
    dst->len = src->len;
}

static void pixie_model_init(PixieModel *m) {
    m->vars = NULL;
    m->n_vars = 0;
    m->cap_vars = 0;
    m->cons = NULL;
    m->n_cons = 0;
    m->cap_cons = 0;
    m->obj_sense = +1;
}

static void pixie_model_free(PixieModel *m) {
    int i;
    for (i = 0; i < m->n_vars; ++i) {
        free(m->vars[i].name);
    }
    for (i = 0; i < m->n_cons; ++i) {
        pixie_sparse_free(&m->cons[i].a);
    }
    free(m->vars);
    free(m->cons);
    m->vars = NULL;
    m->cons = NULL;
    m->n_vars = 0;
    m->cap_vars = 0;
    m->n_cons = 0;
    m->cap_cons = 0;
    m->obj_sense = +1;
}

static int pixie_model_find_var(const PixieModel *m, const char *name) {
    int i;
    for (i = 0; i < m->n_vars; ++i) {
        if (strcmp(m->vars[i].name, name) == 0) {
            return i;
        }
    }
    return -1;
}

static int pixie_model_add_var(PixieModel *m, const char *name) {
    int nc;
    PixieVar *v;
    if (m->n_vars >= m->cap_vars) {
        nc = (m->cap_vars <= 0) ? 8 : m->cap_vars;
        while (nc <= m->n_vars) {
            if (nc > INT_MAX / 2) {
                nc = m->n_vars + 1;
                break;
            }
            nc *= 2;
        }
        m->vars = (PixieVar *)pixie_xrealloc(m->vars, (size_t)nc * sizeof(PixieVar));
        m->cap_vars = nc;
    }
    v = &m->vars[m->n_vars];
    v->name = pixie_strdup_c(name);
    v->obj = 0.0;
    v->lb = 0.0;
    v->ub = HUGE_VAL;
    v->type = PIXIE_VAR_CONT;
    ++m->n_vars;
    return m->n_vars - 1;
}

static int pixie_model_get_var(PixieModel *m, const char *name, bool create) {
    int idx = pixie_model_find_var(m, name);
    if (idx >= 0) {
        return idx;
    }
    if (!create) {
        return -1;
    }
    return pixie_model_add_var(m, name);
}

static int pixie_model_reserve_cons(PixieModel *m, int need) {
    int nc;
    if (need <= m->cap_cons) {
        return 1;
    }
    nc = (m->cap_cons <= 0) ? 8 : m->cap_cons;
    while (nc < need) {
        if (nc > INT_MAX / 2) {
            nc = need;
            break;
        }
        nc *= 2;
    }
    m->cons = (PixieConstraint *)pixie_xrealloc(m->cons, (size_t)nc * sizeof(PixieConstraint));
    m->cap_cons = nc;
    return 1;
}

static int pixie_model_add_empty_constraint(PixieModel *m, PixieCompare cmp, double rhs) {
    if (!pixie_model_reserve_cons(m, m->n_cons + 1)) {
        return -1;
    }
    pixie_sparse_init(&m->cons[m->n_cons].a);
    m->cons[m->n_cons].cmp = cmp;
    m->cons[m->n_cons].rhs = rhs;
    ++m->n_cons;
    return m->n_cons - 1;
}

static int pixie_model_add_constraint(PixieModel *m, const PixieSparse *a, PixieCompare cmp, double rhs) {
    int idx = pixie_model_add_empty_constraint(m, cmp, rhs);
    if (idx < 0) {
        return 0;
    }
    pixie_sparse_clone(&m->cons[idx].a, a);
    pixie_sparse_prune(&m->cons[idx].a, PIXIE_EPS);
    return 1;
}

static void pixie_var_mark_integer(PixieModel *m, int idx) {
    if (idx < 0 || idx >= m->n_vars) {
        return;
    }
    if (m->vars[idx].type == PIXIE_VAR_CONT) {
        m->vars[idx].type = PIXIE_VAR_INT;
    }
}

static void pixie_var_mark_binary(PixieModel *m, int idx) {
    if (idx < 0 || idx >= m->n_vars) {
        return;
    }
    m->vars[idx].type = PIXIE_VAR_BIN;
    if (m->vars[idx].lb < 0.0) {
        m->vars[idx].lb = 0.0;
    }
    if (m->vars[idx].ub > 1.0) {
        m->vars[idx].ub = 1.0;
    }
}

static int pixie_count_integer_vars(const PixieModel *m) {
    int i;
    int cnt = 0;
    for (i = 0; i < m->n_vars; ++i) {
        if (m->vars[i].type != PIXIE_VAR_CONT) {
            ++cnt;
        }
    }
    return cnt;
}

static bool pixie_is_name_start_char(char c) {
    return isalpha((unsigned char)c) || c == '_' || c == '.' || c == '$';
}

static bool pixie_is_name_char(char c) {
    return isalnum((unsigned char)c) || c == '_' || c == '.' || c == '$' || c == '[' || c == ']' || c == '#';
}

static int pixie_find_comparator(const char *s, int *pos, int *len, PixieCompare *cmp) {
    int i;
    int n = (int)strlen(s);
    for (i = 0; i < n; ++i) {
        if (s[i] == '<') {
            if (i + 1 < n && s[i + 1] == '=') {
                *pos = i;
                *len = 2;
                *cmp = PIXIE_CMP_LE;
                return 1;
            }
            *pos = i;
            *len = 1;
            *cmp = PIXIE_CMP_LE;
            return 1;
        }
        if (s[i] == '>') {
            if (i + 1 < n && s[i + 1] == '=') {
                *pos = i;
                *len = 2;
                *cmp = PIXIE_CMP_GE;
                return 1;
            }
            *pos = i;
            *len = 1;
            *cmp = PIXIE_CMP_GE;
            return 1;
        }
        if (s[i] == '=') {
            *pos = i;
            *len = 1;
            *cmp = PIXIE_CMP_EQ;
            return 1;
        }
    }
    return 0;
}

static int pixie_parse_linear_expr(PixieModel *m, const char *text, bool create_vars, PixieSparse *out_terms, double *out_const, char *err, size_t errsz) {
    const char *p = text;
    int sign = +1;
    pixie_sparse_init(out_terms);
    *out_const = 0.0;

    while (*p != '\0') {
        char *end_num = NULL;
        double coef = 0.0;
        bool has_num = false;

        while (*p != '\0' && isspace((unsigned char)*p)) {
            ++p;
        }
        if (*p == '\0') {
            break;
        }
        if (*p == '+') {
            sign = +1;
            ++p;
            continue;
        }
        if (*p == '-') {
            sign = -1;
            ++p;
            continue;
        }

        errno = 0;
        coef = strtod(p, &end_num);
        if (end_num != p) {
            has_num = true;
            if (errno == ERANGE) {
                pixie_sparse_free(out_terms);
                pixie_set_error(err, errsz, "numeric overflow while parsing expression");
                return 0;
            }
            p = end_num;
            while (*p != '\0' && isspace((unsigned char)*p)) {
                ++p;
            }
            if (*p == '*') {
                ++p;
                while (*p != '\0' && isspace((unsigned char)*p)) {
                    ++p;
                }
            }
        }

        if (pixie_is_name_start_char(*p)) {
            const char *start = p;
            char name_buf[256];
            size_t name_len = 0;
            int vidx;

            while (*p != '\0' && pixie_is_name_char(*p)) {
                ++p;
            }
            name_len = (size_t)(p - start);
            if (name_len == 0 || name_len >= sizeof(name_buf)) {
                pixie_sparse_free(out_terms);
                pixie_set_error(err, errsz, "invalid variable name in expression");
                return 0;
            }
            memcpy(name_buf, start, name_len);
            name_buf[name_len] = '\0';
            vidx = pixie_model_get_var(m, name_buf, create_vars);
            if (vidx < 0) {
                pixie_sparse_free(out_terms);
                pixie_set_error(err, errsz, "unknown variable in expression");
                return 0;
            }
            if (!has_num) {
                coef = 1.0;
            }
            pixie_sparse_add(out_terms, vidx, (double)sign * coef);
            sign = +1;
            continue;
        }

        if (has_num) {
            *out_const += (double)sign * coef;
            sign = +1;
            continue;
        }

        {
            char msg[256];
            snprintf(msg, sizeof(msg), "cannot parse linear expression near '%.80s'", p);
            pixie_sparse_free(out_terms);
            pixie_set_error(err, errsz, msg);
        }
        return 0;
    }

    pixie_sparse_prune(out_terms, PIXIE_EPS);
    return 1;
}

/* Portable fallback for systems without strtok_s. */
#ifndef _WIN32
#define strtok_s strtok_r
#endif

static int pixie_parse_lp_objective(PixieModel *m, char *stmt, char *err, size_t errsz) {
    PixieSparse terms;
    double cst = 0.0;
    char *p = stmt;
    char *colon = NULL;
    int i;

    p = pixie_trim(p);
    colon = strchr(p, ':');
    if (colon != NULL) {
        *colon = '\0';
        p = pixie_trim(colon + 1);
    }

    if (!pixie_parse_linear_expr(m, p, true, &terms, &cst, err, errsz)) {
        return 0;
    }
    (void)cst;
    for (i = 0; i < terms.len; ++i) {
        int v = terms.idx[i];
        m->vars[v].obj += terms.val[i];
    }
    pixie_sparse_free(&terms);
    return 1;
}

static int pixie_parse_lp_constraint(PixieModel *m, char *stmt, char *err, size_t errsz) {
    PixieCompare cmp = PIXIE_CMP_EQ;
    char *buf = pixie_strdup_c(stmt);
    char *lhs = NULL;
    char *rhs = NULL;
    char *p = NULL;
    char *colon = NULL;
    int clen = 0;
    PixieSparse lterms;
    PixieSparse rterms;
    PixieSparse comb;
    double lc = 0.0;
    double rc = 0.0;
    double nrhs = 0.0;
    int i;

    lhs = buf;
    p = buf;
    while (*p != '\0') {
        if (*p == '<') {
            cmp = PIXIE_CMP_LE;
            clen = (p[1] == '=') ? 2 : 1;
            break;
        }
        if (*p == '>') {
            cmp = PIXIE_CMP_GE;
            clen = (p[1] == '=') ? 2 : 1;
            break;
        }
        if (*p == '=') {
            cmp = PIXIE_CMP_EQ;
            clen = 1;
            break;
        }
        ++p;
    }
    if (*p == '\0') {
        free(buf);
        pixie_set_error(err, errsz, "constraint has no comparator");
        return 0;
    }
    *p = '\0';
    rhs = p + clen;

    lhs = pixie_trim(lhs);
    rhs = pixie_trim(rhs);

    colon = strchr(lhs, ':');
    if (colon != NULL) {
        lhs = pixie_trim(colon + 1);
    }

    if (strchr(lhs, '<') != NULL || strchr(lhs, '>') != NULL || strchr(lhs, '=') != NULL) {
        char msg[256];
        snprintf(msg, sizeof(msg), "internal split error lhs='%.120s'", lhs);
        free(buf);
        pixie_set_error(err, errsz, msg);
        return 0;
    }

    if (!pixie_parse_linear_expr(m, lhs, true, &lterms, &lc, err, errsz)) {
        free(buf);
        return 0;
    }
    if (!pixie_parse_linear_expr(m, rhs, true, &rterms, &rc, err, errsz)) {
        pixie_sparse_free(&lterms);
        free(buf);
        return 0;
    }

    pixie_sparse_init(&comb);
    for (i = 0; i < lterms.len; ++i) {
        pixie_sparse_add(&comb, lterms.idx[i], lterms.val[i]);
    }
    for (i = 0; i < rterms.len; ++i) {
        pixie_sparse_add(&comb, rterms.idx[i], -rterms.val[i]);
    }
    nrhs = rc - lc;
    pixie_sparse_prune(&comb, PIXIE_EPS);

    if (!pixie_model_add_constraint(m, &comb, cmp, nrhs)) {
        pixie_sparse_free(&comb);
        pixie_sparse_free(&lterms);
        pixie_sparse_free(&rterms);
        free(buf);
        pixie_set_error(err, errsz, "failed to add LP constraint");
        return 0;
    }

    pixie_sparse_free(&comb);
    pixie_sparse_free(&lterms);
    pixie_sparse_free(&rterms);
    free(buf);
    return 1;
}

static void pixie_tokvec_init(PixieTokVec *tv) {
    tv->data = NULL;
    tv->len = 0;
    tv->cap = 0;
}

static void pixie_tokvec_free(PixieTokVec *tv) {
    int i;
    for (i = 0; i < tv->len; ++i) {
        free(tv->data[i]);
    }
    free(tv->data);
    tv->data = NULL;
    tv->len = 0;
    tv->cap = 0;
}

static void pixie_tokvec_push(PixieTokVec *tv, const char *tok) {
    int nc;
    if (tv->len >= tv->cap) {
        nc = (tv->cap <= 0) ? 8 : tv->cap;
        while (nc <= tv->len) {
            if (nc > INT_MAX / 2) {
                nc = tv->len + 1;
                break;
            }
            nc *= 2;
        }
        tv->data = (char **)pixie_xrealloc(tv->data, (size_t)nc * sizeof(char *));
        tv->cap = nc;
    }
    tv->data[tv->len] = pixie_strdup_c(tok);
    ++tv->len;
}

static void pixie_tokenize_bounds(const char *stmt, PixieTokVec *tv) {
    const char *p = stmt;
    char buf[512];
    pixie_tokvec_init(tv);

    while (*p != '\0') {
        size_t n = 0;
        while (*p != '\0' && isspace((unsigned char)*p)) {
            ++p;
        }
        if (*p == '\0') {
            break;
        }
        if (*p == '<' || *p == '>' || *p == '=') {
            buf[n++] = *p;
            if ((p[0] == '<' || p[0] == '>') && p[1] == '=') {
                ++p;
                buf[n++] = *p;
            }
            ++p;
            buf[n] = '\0';
            pixie_tokvec_push(tv, buf);
            continue;
        }
        while (*p != '\0' && !isspace((unsigned char)*p) && *p != '<' && *p != '>' && *p != '=') {
            if (n + 1 < sizeof(buf)) {
                buf[n++] = *p;
            }
            ++p;
        }
        buf[n] = '\0';
        if (n > 0) {
            pixie_tokvec_push(tv, buf);
        }
    }
}

static int pixie_lp_apply_bound_rel(PixieModel *m, const char *left, const char *op, const char *right, char *err, size_t errsz) {
    double lv = 0.0;
    double rv = 0.0;
    bool left_num = pixie_parse_double_token(left, &lv);
    bool right_num = pixie_parse_double_token(right, &rv);
    int vidx = -1;

    if (!left_num && right_num) {
        vidx = pixie_model_get_var(m, left, true);
        if (vidx < 0) {
            pixie_set_error(err, errsz, "invalid bound variable");
            return 0;
        }
        if (strcmp(op, "<=") == 0 || strcmp(op, "<") == 0) {
            m->vars[vidx].ub = pixie_min(m->vars[vidx].ub, rv);
            return 1;
        }
        if (strcmp(op, ">=") == 0 || strcmp(op, ">") == 0) {
            m->vars[vidx].lb = pixie_max(m->vars[vidx].lb, rv);
            return 1;
        }
        if (strcmp(op, "=") == 0) {
            m->vars[vidx].lb = rv;
            m->vars[vidx].ub = rv;
            return 1;
        }
        pixie_set_error(err, errsz, "unsupported bounds operator");
        return 0;
    }

    if (left_num && !right_num) {
        vidx = pixie_model_get_var(m, right, true);
        if (vidx < 0) {
            pixie_set_error(err, errsz, "invalid bound variable");
            return 0;
        }
        if (strcmp(op, "<=") == 0 || strcmp(op, "<") == 0) {
            m->vars[vidx].lb = pixie_max(m->vars[vidx].lb, lv);
            return 1;
        }
        if (strcmp(op, ">=") == 0 || strcmp(op, ">") == 0) {
            m->vars[vidx].ub = pixie_min(m->vars[vidx].ub, lv);
            return 1;
        }
        if (strcmp(op, "=") == 0) {
            m->vars[vidx].lb = lv;
            m->vars[vidx].ub = lv;
            return 1;
        }
        pixie_set_error(err, errsz, "unsupported bounds operator");
        return 0;
    }

    pixie_set_error(err, errsz, "bounds relation must be between one variable and one number");
    return 0;
}

static int pixie_parse_lp_bounds(PixieModel *m, char *stmt, char *err, size_t errsz) {
    PixieTokVec tv;
    char *t;
    int idx;
    pixie_tokenize_bounds(stmt, &tv);
    if (tv.len == 0) {
        pixie_tokvec_free(&tv);
        return 1;
    }

    t = tv.data[0];
    if (tv.len >= 2 && pixie_stricmp_c(tv.data[1], "free") == 0) {
        idx = pixie_model_get_var(m, t, true);
        if (idx < 0) {
            pixie_tokvec_free(&tv);
            pixie_set_error(err, errsz, "invalid variable in free bound");
            return 0;
        }
        m->vars[idx].lb = -HUGE_VAL;
        m->vars[idx].ub = HUGE_VAL;
        pixie_tokvec_free(&tv);
        return 1;
    }

    if (tv.len == 3) {
        int ok = pixie_lp_apply_bound_rel(m, tv.data[0], tv.data[1], tv.data[2], err, errsz);
        pixie_tokvec_free(&tv);
        return ok;
    }
    if (tv.len == 5) {
        int ok1 = pixie_lp_apply_bound_rel(m, tv.data[0], tv.data[1], tv.data[2], err, errsz);
        int ok2 = 0;
        if (!ok1) {
            pixie_tokvec_free(&tv);
            return 0;
        }
        ok2 = pixie_lp_apply_bound_rel(m, tv.data[2], tv.data[3], tv.data[4], err, errsz);
        pixie_tokvec_free(&tv);
        return ok2;
    }

    pixie_tokvec_free(&tv);
    pixie_set_error(err, errsz, "unsupported bounds syntax");
    return 0;
}

static bool pixie_lp_is_keyword(const char *s) {
    return pixie_stricmp_c(s, "int") == 0 ||
           pixie_stricmp_c(s, "integer") == 0 ||
           pixie_stricmp_c(s, "integers") == 0 ||
           pixie_stricmp_c(s, "general") == 0 ||
           pixie_stricmp_c(s, "generals") == 0 ||
           pixie_stricmp_c(s, "binary") == 0 ||
           pixie_stricmp_c(s, "binaries") == 0 ||
           pixie_stricmp_c(s, "bin") == 0;
}

static int pixie_parse_lp_var_list(PixieModel *m, char *stmt, bool binary, char *err, size_t errsz) {
    char *dup = pixie_strdup_c(stmt);
    char *tok = NULL;
    char *ctx = NULL;
    tok = strtok_s(dup, " \t\r\n,:", &ctx);
    while (tok != NULL) {
        int idx;
        if (!pixie_lp_is_keyword(tok)) {
            idx = pixie_model_get_var(m, tok, true);
            if (idx < 0) {
                free(dup);
                pixie_set_error(err, errsz, "invalid variable in integer/binary list");
                return 0;
            }
            if (binary) {
                pixie_var_mark_binary(m, idx);
            } else {
                pixie_var_mark_integer(m, idx);
            }
        }
        tok = strtok_s(NULL, " \t\r\n,:", &ctx);
    }
    free(dup);
    return 1;
}

enum {
    LP_SEC_NONE = 0,
    LP_SEC_OBJ = 1,
    LP_SEC_CONS = 2,
    LP_SEC_BOUNDS = 3,
    LP_SEC_GENERAL = 4,
    LP_SEC_BINARY = 5
};

static int pixie_parse_lp_statement(PixieModel *m, char *stmt, int *section, char *err, size_t errsz) {
    char *s = pixie_trim(stmt);
    int pos = -1;
    int clen = 0;
    PixieCompare cmp = PIXIE_CMP_EQ;

    if (*s == '\0') {
        return 1;
    }

    if (pixie_starts_with_ci(s, "min:")) {
        m->obj_sense = +1;
        *section = LP_SEC_OBJ;
        s = pixie_trim(s + 4);
        if (*s == '\0') {
            return 1;
        }
        return pixie_parse_lp_objective(m, s, err, errsz);
    }
    if (pixie_starts_with_ci(s, "max:")) {
        m->obj_sense = -1;
        *section = LP_SEC_OBJ;
        s = pixie_trim(s + 4);
        if (*s == '\0') {
            return 1;
        }
        return pixie_parse_lp_objective(m, s, err, errsz);
    }
    if (pixie_starts_with_ci(s, "minimize")) {
        m->obj_sense = +1;
        *section = LP_SEC_OBJ;
        s = pixie_trim(s + 8);
        if (*s == '\0') {
            return 1;
        }
        return pixie_parse_lp_objective(m, s, err, errsz);
    }
    if (pixie_starts_with_ci(s, "maximize")) {
        m->obj_sense = -1;
        *section = LP_SEC_OBJ;
        s = pixie_trim(s + 8);
        if (*s == '\0') {
            return 1;
        }
        return pixie_parse_lp_objective(m, s, err, errsz);
    }
    if (pixie_starts_with_ci(s, "subject to")) {
        *section = LP_SEC_CONS;
        s = pixie_trim(s + 10);
        if (*s == '\0') {
            return 1;
        }
    }
    if (pixie_starts_with_ci(s, "such that")) {
        *section = LP_SEC_CONS;
        s = pixie_trim(s + 9);
        if (*s == '\0') {
            return 1;
        }
    }
    if (pixie_starts_with_ci(s, "st")) {
        if (s[2] == '\0' || isspace((unsigned char)s[2])) {
            *section = LP_SEC_CONS;
            s = pixie_trim(s + 2);
            if (*s == '\0') {
                return 1;
            }
        }
    }
    if (pixie_stricmp_c(s, "bounds") == 0 || pixie_starts_with_ci(s, "bounds ")) {
        *section = LP_SEC_BOUNDS;
        s = pixie_trim(s + 6);
        if (*s == '\0') {
            return 1;
        }
        return pixie_parse_lp_bounds(m, s, err, errsz);
    }
    if (pixie_stricmp_c(s, "general") == 0 || pixie_stricmp_c(s, "generals") == 0 ||
        pixie_stricmp_c(s, "integer") == 0 || pixie_stricmp_c(s, "integers") == 0 ||
        pixie_stricmp_c(s, "int") == 0 || pixie_starts_with_ci(s, "general ") ||
        pixie_starts_with_ci(s, "generals ") || pixie_starts_with_ci(s, "integer ") ||
        pixie_starts_with_ci(s, "integers ") || pixie_starts_with_ci(s, "int ")) {
        *section = LP_SEC_GENERAL;
        return pixie_parse_lp_var_list(m, s, false, err, errsz);
    }
    if (pixie_stricmp_c(s, "binary") == 0 || pixie_stricmp_c(s, "binaries") == 0 ||
        pixie_stricmp_c(s, "bin") == 0 || pixie_starts_with_ci(s, "binary ") ||
        pixie_starts_with_ci(s, "binaries ") || pixie_starts_with_ci(s, "bin ")) {
        *section = LP_SEC_BINARY;
        return pixie_parse_lp_var_list(m, s, true, err, errsz);
    }
    if (pixie_stricmp_c(s, "end") == 0) {
        return 2;
    }

    if (pixie_find_comparator(s, &pos, &clen, &cmp)) {
        *section = LP_SEC_CONS;
        return pixie_parse_lp_constraint(m, s, err, errsz);
    }

    if (*section == LP_SEC_OBJ) {
        return pixie_parse_lp_objective(m, s, err, errsz);
    }
    if (*section == LP_SEC_BOUNDS) {
        return pixie_parse_lp_bounds(m, s, err, errsz);
    }
    if (*section == LP_SEC_GENERAL) {
        return pixie_parse_lp_var_list(m, s, false, err, errsz);
    }
    if (*section == LP_SEC_BINARY) {
        return pixie_parse_lp_var_list(m, s, true, err, errsz);
    }

    if (*section == LP_SEC_CONS) {
        pixie_set_error(err, errsz, "constraint line missing comparator");
        return 0;
    }

    return pixie_parse_lp_constraint(m, s, err, errsz);
}

static int pixie_parse_lp_file(const char *path, PixieModel *m, char *err, size_t errsz) {
    FILE *fp = fopen(path, "rb");
    char line[4096];
    int section = LP_SEC_NONE;
    int ended = 0;
    if (fp == NULL) {
        pixie_set_error(err, errsz, "failed to open LP file");
        return 0;
    }

    while (fgets(line, (int)sizeof(line), fp) != NULL) {
        char *cursor = line;
        char *comment = strchr(cursor, '\\');
        int line_done = 0;
        if (comment != NULL) {
            *comment = '\0';
        }
        while (!line_done) {
            char *semi = strchr(cursor, ';');
            char *stmt = NULL;
            int rc;
            if (semi != NULL) {
                *semi = '\0';
                stmt = pixie_trim(cursor);
                if (*stmt != '\0') {
                    rc = pixie_parse_lp_statement(m, stmt, &section, err, errsz);
                    if (rc == 0) {
                        if (err != NULL && errsz > 0) {
                            char msg[512];
                            snprintf(msg, sizeof(msg), "%s | %.240s", err, stmt);
                            pixie_set_error(err, errsz, msg);
                        }
                        fclose(fp);
                        return 0;
                    }
                    if (rc == 2) {
                        ended = 1;
                        break;
                    }
                }
                cursor = semi + 1;
            } else {
                stmt = pixie_trim(cursor);
                if (*stmt != '\0') {
                    rc = pixie_parse_lp_statement(m, stmt, &section, err, errsz);
                    if (rc == 0) {
                        if (err != NULL && errsz > 0) {
                            char msg[512];
                            snprintf(msg, sizeof(msg), "%s | %.240s", err, stmt);
                            pixie_set_error(err, errsz, msg);
                        }
                        fclose(fp);
                        return 0;
                    }
                    if (rc == 2) {
                        ended = 1;
                        break;
                    }
                }
                line_done = 1;
            }
        }
        if (ended) {
            break;
        }
    }

    fclose(fp);
    if (m->n_vars <= 0) {
        pixie_set_error(err, errsz, "LP parser found no variables");
        return 0;
    }
    return 1;
}

static int pixie_split_ws(char *line, char **tok, int max_tok) {
    int n = 0;
    char *p = line;
    while (*p != '\0') {
        while (*p != '\0' && isspace((unsigned char)*p)) {
            ++p;
        }
        if (*p == '\0') {
            break;
        }
        if (n >= max_tok) {
            break;
        }
        tok[n++] = p;
        while (*p != '\0' && !isspace((unsigned char)*p)) {
            ++p;
        }
        if (*p == '\0') {
            break;
        }
        *p = '\0';
        ++p;
    }
    return n;
}

typedef struct {
    char *name;
    int cidx;
    PixieCompare cmp;
} PixieMpsRow;

typedef struct {
    PixieMpsRow *data;
    int len;
    int cap;
} PixieMpsRows;

static void pixie_mps_rows_init(PixieMpsRows *r) {
    r->data = NULL;
    r->len = 0;
    r->cap = 0;
}

static void pixie_mps_rows_free(PixieMpsRows *r) {
    int i;
    for (i = 0; i < r->len; ++i) {
        free(r->data[i].name);
    }
    free(r->data);
    r->data = NULL;
    r->len = 0;
    r->cap = 0;
}

static int pixie_mps_rows_add(PixieMpsRows *r, const char *name, int cidx, PixieCompare cmp) {
    int nc;
    if (r->len >= r->cap) {
        nc = (r->cap <= 0) ? 8 : r->cap;
        while (nc <= r->len) {
            if (nc > INT_MAX / 2) {
                nc = r->len + 1;
                break;
            }
            nc *= 2;
        }
        r->data = (PixieMpsRow *)pixie_xrealloc(r->data, (size_t)nc * sizeof(PixieMpsRow));
        r->cap = nc;
    }
    r->data[r->len].name = pixie_strdup_c(name);
    r->data[r->len].cidx = cidx;
    r->data[r->len].cmp = cmp;
    ++r->len;
    return r->len - 1;
}

static int pixie_mps_rows_find(const PixieMpsRows *r, const char *name) {
    int i;
    for (i = 0; i < r->len; ++i) {
        if (strcmp(r->data[i].name, name) == 0) {
            return i;
        }
    }
    return -1;
}

static int pixie_parse_mps_file(const char *path, PixieModel *m, char *err, size_t errsz) {
    FILE *fp = fopen(path, "rb");
    char line[4096];
    char *tok[32];
    int ntok;
    enum { MPS_NONE = 0, MPS_ROWS = 1, MPS_COLUMNS = 2, MPS_RHS = 3, MPS_BOUNDS = 4, MPS_OBJSENSE = 5 } sec = MPS_NONE;
    PixieMpsRows rows;
    int obj_row = -1;
    bool int_mode = false;
    int i;

    if (fp == NULL) {
        pixie_set_error(err, errsz, "failed to open MPS file");
        return 0;
    }

    pixie_mps_rows_init(&rows);
    m->obj_sense = +1;

    while (fgets(line, (int)sizeof(line), fp) != NULL) {
        bool indented = (line[0] != '\0' && isspace((unsigned char)line[0])) ? true : false;
        char *s = pixie_trim(line);
        if (*s == '\0') {
            continue;
        }
        if (*s == '*') {
            if (strstr(s, "origsense='MAX'") != NULL || strstr(s, "origsense=\"MAX\"") != NULL) {
                m->obj_sense = -1;
            } else if (strstr(s, "origsense='MIN'") != NULL || strstr(s, "origsense=\"MIN\"") != NULL) {
                m->obj_sense = +1;
            }
            continue;
        }

        ntok = pixie_split_ws(s, tok, 32);
        if (ntok <= 0) {
            continue;
        }

        if (!indented) {
            if (pixie_stricmp_c(tok[0], "NAME") == 0) {
                sec = MPS_NONE;
                continue;
            }
            if (pixie_stricmp_c(tok[0], "OBJSENSE") == 0) {
                sec = MPS_OBJSENSE;
                continue;
            }
            if (pixie_stricmp_c(tok[0], "ROWS") == 0) {
                sec = MPS_ROWS;
                continue;
            }
            if (pixie_stricmp_c(tok[0], "COLUMNS") == 0) {
                sec = MPS_COLUMNS;
                continue;
            }
            if (pixie_stricmp_c(tok[0], "RHS") == 0) {
                sec = MPS_RHS;
                continue;
            }
            if (pixie_stricmp_c(tok[0], "BOUNDS") == 0) {
                sec = MPS_BOUNDS;
                continue;
            }
            if (pixie_stricmp_c(tok[0], "ENDATA") == 0 || pixie_stricmp_c(tok[0], "END") == 0) {
                break;
            }
        }

        if (sec == MPS_OBJSENSE) {
            if (pixie_stricmp_c(tok[0], "MAX") == 0 || pixie_stricmp_c(tok[0], "MAXIMIZE") == 0) {
                m->obj_sense = -1;
            } else if (pixie_stricmp_c(tok[0], "MIN") == 0 || pixie_stricmp_c(tok[0], "MINIMIZE") == 0) {
                m->obj_sense = +1;
            }
            continue;
        }

        if (sec == MPS_ROWS) {
            char t;
            int cidx;
            if (ntok < 2) {
                pixie_set_error(err, errsz, "invalid ROWS line in MPS");
                pixie_mps_rows_free(&rows);
                fclose(fp);
                return 0;
            }
            t = tok[0][0];
            if (t == 'N' || t == 'n') {
                if (obj_row < 0) {
                    obj_row = pixie_mps_rows_add(&rows, tok[1], -1, PIXIE_CMP_EQ);
                } else {
                    pixie_mps_rows_add(&rows, tok[1], -1, PIXIE_CMP_EQ);
                }
                continue;
            }
            if (t == 'L' || t == 'l') {
                cidx = pixie_model_add_empty_constraint(m, PIXIE_CMP_LE, 0.0);
            } else if (t == 'G' || t == 'g') {
                cidx = pixie_model_add_empty_constraint(m, PIXIE_CMP_GE, 0.0);
            } else if (t == 'E' || t == 'e') {
                cidx = pixie_model_add_empty_constraint(m, PIXIE_CMP_EQ, 0.0);
            } else {
                pixie_set_error(err, errsz, "unsupported ROWS type in MPS");
                pixie_mps_rows_free(&rows);
                fclose(fp);
                return 0;
            }
            pixie_mps_rows_add(&rows, tok[1], cidx, m->cons[cidx].cmp);
            continue;
        }

        if (sec == MPS_COLUMNS) {
            int vidx;
            if (ntok >= 3 && strcmp(tok[1], "'MARKER'") == 0) {
                if (ntok >= 3 && (strstr(tok[2], "INTORG") != NULL || strstr(tok[2], "intorg") != NULL)) {
                    int_mode = true;
                } else if (ntok >= 3 && (strstr(tok[2], "INTEND") != NULL || strstr(tok[2], "intend") != NULL)) {
                    int_mode = false;
                }
                continue;
            }
            if (ntok < 3) {
                pixie_set_error(err, errsz, "invalid COLUMNS line in MPS");
                pixie_mps_rows_free(&rows);
                fclose(fp);
                return 0;
            }
            vidx = pixie_model_get_var(m, tok[0], true);
            if (vidx < 0) {
                pixie_set_error(err, errsz, "failed to create variable in MPS");
                pixie_mps_rows_free(&rows);
                fclose(fp);
                return 0;
            }
            if (int_mode) {
                pixie_var_mark_integer(m, vidx);
            }
            for (i = 1; i + 1 < ntok; i += 2) {
                int ridx = pixie_mps_rows_find(&rows, tok[i]);
                double v = 0.0;
                if (ridx < 0) {
                    pixie_set_error(err, errsz, "unknown row name in COLUMNS");
                    pixie_mps_rows_free(&rows);
                    fclose(fp);
                    return 0;
                }
                if (!pixie_parse_double_token(tok[i + 1], &v)) {
                    pixie_set_error(err, errsz, "invalid numeric value in COLUMNS");
                    pixie_mps_rows_free(&rows);
                    fclose(fp);
                    return 0;
                }
                if (rows.data[ridx].cidx < 0) {
                    m->vars[vidx].obj += v;
                } else {
                    pixie_sparse_add(&m->cons[rows.data[ridx].cidx].a, vidx, v);
                }
            }
            continue;
        }

        if (sec == MPS_RHS) {
            if (ntok < 3) {
                pixie_set_error(err, errsz, "invalid RHS line in MPS");
                pixie_mps_rows_free(&rows);
                fclose(fp);
                return 0;
            }
            for (i = 1; i + 1 < ntok; i += 2) {
                int ridx = pixie_mps_rows_find(&rows, tok[i]);
                double v = 0.0;
                if (ridx < 0) {
                    pixie_set_error(err, errsz, "unknown row in RHS");
                    pixie_mps_rows_free(&rows);
                    fclose(fp);
                    return 0;
                }
                if (!pixie_parse_double_token(tok[i + 1], &v)) {
                    pixie_set_error(err, errsz, "invalid RHS value");
                    pixie_mps_rows_free(&rows);
                    fclose(fp);
                    return 0;
                }
                if (rows.data[ridx].cidx >= 0) {
                    m->cons[rows.data[ridx].cidx].rhs = v;
                }
            }
            continue;
        }

        if (sec == MPS_BOUNDS) {
            int vidx;
            double v = 0.0;
            if (ntok < 3) {
                pixie_set_error(err, errsz, "invalid BOUNDS line in MPS");
                pixie_mps_rows_free(&rows);
                fclose(fp);
                return 0;
            }
            vidx = pixie_model_get_var(m, tok[2], true);
            if (vidx < 0) {
                pixie_set_error(err, errsz, "unknown variable in BOUNDS");
                pixie_mps_rows_free(&rows);
                fclose(fp);
                return 0;
            }
            if (ntok >= 4) {
                if (!pixie_parse_double_token(tok[3], &v)) {
                    pixie_set_error(err, errsz, "invalid numeric value in BOUNDS");
                    pixie_mps_rows_free(&rows);
                    fclose(fp);
                    return 0;
                }
            }
            if (pixie_stricmp_c(tok[0], "LO") == 0) {
                m->vars[vidx].lb = pixie_max(m->vars[vidx].lb, v);
            } else if (pixie_stricmp_c(tok[0], "UP") == 0) {
                m->vars[vidx].ub = pixie_min(m->vars[vidx].ub, v);
            } else if (pixie_stricmp_c(tok[0], "FX") == 0) {
                m->vars[vidx].lb = v;
                m->vars[vidx].ub = v;
            } else if (pixie_stricmp_c(tok[0], "FR") == 0) {
                m->vars[vidx].lb = -HUGE_VAL;
                m->vars[vidx].ub = HUGE_VAL;
            } else if (pixie_stricmp_c(tok[0], "MI") == 0) {
                m->vars[vidx].lb = -HUGE_VAL;
            } else if (pixie_stricmp_c(tok[0], "PL") == 0) {
                m->vars[vidx].ub = HUGE_VAL;
            } else if (pixie_stricmp_c(tok[0], "BV") == 0) {
                pixie_var_mark_binary(m, vidx);
            } else if (pixie_stricmp_c(tok[0], "LI") == 0) {
                pixie_var_mark_integer(m, vidx);
                m->vars[vidx].lb = pixie_max(m->vars[vidx].lb, v);
            } else if (pixie_stricmp_c(tok[0], "UI") == 0) {
                pixie_var_mark_integer(m, vidx);
                m->vars[vidx].ub = pixie_min(m->vars[vidx].ub, v);
            } else {
                pixie_set_error(err, errsz, "unsupported BOUNDS type in MPS");
                pixie_mps_rows_free(&rows);
                fclose(fp);
                return 0;
            }
            continue;
        }
    }

    fclose(fp);
    for (i = 0; i < m->n_cons; ++i) {
        pixie_sparse_prune(&m->cons[i].a, PIXIE_EPS);
    }
    pixie_mps_rows_free(&rows);
    if (m->n_vars <= 0) {
        pixie_set_error(err, errsz, "MPS parser found no variables");
        return 0;
    }
    return 1;
}

static int pixie_read_model(const char *path, PixieFormat fmt, PixieModel *m, char *err, size_t errsz) {
    int ok = 0;
    if (fmt == PIXIE_FORMAT_LP) {
        return pixie_parse_lp_file(path, m, err, errsz);
    }
    if (fmt == PIXIE_FORMAT_MPS) {
        return pixie_parse_mps_file(path, m, err, errsz);
    }
    if (pixie_ends_with_ci(path, ".lp")) {
        return pixie_parse_lp_file(path, m, err, errsz);
    }
    if (pixie_ends_with_ci(path, ".mps")) {
        return pixie_parse_mps_file(path, m, err, errsz);
    }
    ok = pixie_parse_lp_file(path, m, err, errsz);
    if (ok) {
        return 1;
    }
    pixie_model_free(m);
    pixie_model_init(m);
    ok = pixie_parse_mps_file(path, m, err, errsz);
    if (ok) {
        return 1;
    }
    pixie_set_error(err, errsz, "could not parse file as LP or MPS");
    return 0;
}

static void pixie_convec_init(PixieConVec *v) {
    v->data = NULL;
    v->len = 0;
    v->cap = 0;
}

static void pixie_convec_free(PixieConVec *v) {
    int i;
    for (i = 0; i < v->len; ++i) {
        pixie_sparse_free(&v->data[i].a);
    }
    free(v->data);
    v->data = NULL;
    v->len = 0;
    v->cap = 0;
}

static int pixie_convec_push(PixieConVec *v, PixieConstraint *c) {
    int nc;
    if (v->len >= v->cap) {
        nc = (v->cap <= 0) ? 8 : v->cap;
        while (nc <= v->len) {
            if (nc > INT_MAX / 2) {
                nc = v->len + 1;
                break;
            }
            nc *= 2;
        }
        v->data = (PixieConstraint *)pixie_xrealloc(v->data, (size_t)nc * sizeof(PixieConstraint));
        v->cap = nc;
    }
    v->data[v->len].cmp = c->cmp;
    v->data[v->len].rhs = c->rhs;
    pixie_sparse_move(&v->data[v->len].a, &c->a);
    ++v->len;
    return 1;
}

static void pixie_simplex_free(PixieSimplex *sm) {
    free(sm->raw);
    free(sm->table);
    free(sm->v_idx);
    free(sm->a_flg);
    sm->raw = NULL;
    sm->table = NULL;
    sm->v_idx = NULL;
    sm->a_flg = NULL;
    sm->v_cnt = 0;
    sm->e_cnt = 0;
    sm->s_cnt = 0;
    sm->a_cnt = 0;
}

static void pixie_standardize_row(const PixieConstraint *c, bool *flip_sign, bool *has_slack, double *slack_coef, bool *has_artificial) {
    PixieCompare cmp = c->cmp;
    double rhs = c->rhs;
    bool flip = false;

    if (rhs < -PIXIE_FEAS_TOL || (cmp == PIXIE_CMP_GE && rhs <= PIXIE_FEAS_TOL)) {
        flip = true;
        if (cmp == PIXIE_CMP_LE) {
            cmp = PIXIE_CMP_GE;
        } else if (cmp == PIXIE_CMP_GE) {
            cmp = PIXIE_CMP_LE;
        }
    }

    *flip_sign = flip;
    if (cmp == PIXIE_CMP_EQ) {
        *has_slack = false;
        *slack_coef = 0.0;
        *has_artificial = true;
    } else if (cmp == PIXIE_CMP_LE) {
        *has_slack = true;
        *slack_coef = 1.0;
        *has_artificial = false;
    } else {
        *has_slack = true;
        *slack_coef = -1.0;
        *has_artificial = true;
    }
}

static int pixie_count_slack_artificial(const PixieConstraint *cons, int e_cnt, int *s_out, int *a_out) {
    int i;
    int s_cnt = 0;
    int a_cnt = 0;
    for (i = 0; i < e_cnt; ++i) {
        bool flip = false;
        bool has_slack = false;
        bool has_art = false;
        double slack_coef = 0.0;
        pixie_standardize_row(&cons[i], &flip, &has_slack, &slack_coef, &has_art);
        (void)flip;
        (void)slack_coef;
        if (has_slack) {
            ++s_cnt;
        }
        if (has_art) {
            ++a_cnt;
        }
    }
    *s_out = s_cnt;
    *a_out = a_cnt;
    return 1;
}

static int pixie_make_simplex_table(const PixieConstraint *cons, int e_cnt, int v_cnt, const double *obj_max, PixieSimplex *sm) {
    int i;
    int j;
    int k;
    int s_cnt = 0;
    int a_cnt = 0;
    int rows;
    int cols;
    int idx;

    pixie_count_slack_artificial(cons, e_cnt, &s_cnt, &a_cnt);

    rows = e_cnt + 1;
    cols = v_cnt + s_cnt + a_cnt + 1;

    sm->raw = (double *)pixie_xcalloc((size_t)rows * (size_t)cols, sizeof(double));
    sm->table = (double **)pixie_xmalloc((size_t)rows * sizeof(double *));
    for (i = 0; i < rows; ++i) {
        sm->table[i] = sm->raw + (size_t)i * (size_t)cols;
    }
    sm->v_idx = (int *)pixie_xcalloc((size_t)e_cnt, sizeof(int));
    sm->a_flg = (bool *)pixie_xcalloc((size_t)e_cnt, sizeof(bool));
    sm->v_cnt = v_cnt;
    sm->e_cnt = e_cnt;
    sm->s_cnt = s_cnt;
    sm->a_cnt = a_cnt;

    idx = v_cnt + 1;
    for (i = 0; i < e_cnt; ++i) {
        bool flip = false;
        bool has_slack = false;
        bool has_art = false;
        double slack_coef = 0.0;
        pixie_standardize_row(&cons[i], &flip, &has_slack, &slack_coef, &has_art);

        sm->table[i][0] = cons[i].rhs;
        for (k = 0; k < cons[i].a.len; ++k) {
            int v = cons[i].a.idx[k];
            if (v >= 0 && v < v_cnt) {
                sm->table[i][v + 1] += cons[i].a.val[k];
            }
        }
        if (flip) {
            for (j = 0; j <= v_cnt; ++j) {
                sm->table[i][j] = -sm->table[i][j];
            }
        }
        if (has_slack) {
            sm->table[i][idx] = slack_coef;
            if (!has_art && slack_coef > 0.0) {
                sm->v_idx[i] = idx - 1;
            }
            ++idx;
        }
        sm->a_flg[i] = has_art;
    }

    for (i = 0; i < e_cnt; ++i) {
        if (sm->a_flg[i]) {
            sm->table[i][idx] = 1.0;
            sm->v_idx[i] = idx - 1;
            ++idx;
        }
    }

    if (a_cnt == 0) {
        for (i = 0; i < v_cnt; ++i) {
            sm->table[e_cnt][i + 1] = -obj_max[i];
        }
    } else {
        for (i = 0; i < e_cnt; ++i) {
            if (!sm->a_flg[i]) {
                continue;
            }
            for (j = 0; j <= v_cnt + s_cnt; ++j) {
                sm->table[e_cnt][j] -= sm->table[i][j];
            }
        }
    }

    return 1;
}

static void pixie_change_z_param(PixieSimplex *sm, int v_cnt, const double *obj_max) {
    int i;
    int j;
    int s_cnt = sm->s_cnt;
    int e_cnt = sm->e_cnt;
    for (j = 0; j <= v_cnt + s_cnt; ++j) {
        sm->table[e_cnt][j] = 0.0;
        if (j >= 1 && j <= v_cnt) {
            sm->table[e_cnt][j] = -obj_max[j - 1];
        }
        for (i = 0; i < e_cnt; ++i) {
            if (sm->v_idx[i] < v_cnt) {
                sm->table[e_cnt][j] += obj_max[sm->v_idx[i]] * sm->table[i][j];
            } else if (sm->v_idx[i] < v_cnt + s_cnt) {
            }
        }
    }
}

static PixieSimplexStatus pixie_simplex_solve(PixieSimplex *sm, int enter_columns) {
    int simplex_column = sm->v_cnt + sm->s_cnt + sm->a_cnt;
    int total_columns = simplex_column;
    long long iter = 0;
    long long max_iter = 1000000LL + 1000LL * (long long)(sm->e_cnt + 1) * (long long)(simplex_column + 1);
    if (enter_columns <= 0 || enter_columns > total_columns) {
        enter_columns = total_columns;
    }
    while (1) {
        int col = -1;
        int i;
        int j;
        int pivot_row = -1;
        double best_ratio = 0.0;
        int best_basis = INT_MAX;

        if (++iter > max_iter) {
            return PIXIE_SIMPLEX_NUMERIC;
        }

        for (j = 1; j <= enter_columns; ++j) {
            if (sm->table[sm->e_cnt][j] < -PIXIE_PIV_TOL) {
                col = j;
                break;
            }
        }
        if (col < 0) {
            return PIXIE_SIMPLEX_OPTIMAL;
        }

        for (i = 0; i < sm->e_cnt; ++i) {
            double a = sm->table[i][col];
            double rhs = sm->table[i][0];
            double ratio;
            if (a <= PIXIE_PIV_TOL) {
                continue;
            }
            ratio = rhs / a;
            if (pivot_row < 0 ||
                ratio < best_ratio - PIXIE_FEAS_TOL ||
                (fabs(ratio - best_ratio) <= PIXIE_FEAS_TOL && sm->v_idx[i] < best_basis)) {
                pivot_row = i;
                best_ratio = ratio;
                best_basis = sm->v_idx[i];
            }
        }

        if (pivot_row < 0) {
            return PIXIE_SIMPLEX_UNBOUNDED;
        }

        {
            double piv = sm->table[pivot_row][col];
            sm->v_idx[pivot_row] = col - 1;
            for (j = 0; j <= total_columns; ++j) {
                sm->table[pivot_row][j] /= piv;
            }
        }
        for (i = 0; i <= sm->e_cnt; ++i) {
            double factor;
            if (i == pivot_row) {
                continue;
            }
            factor = sm->table[i][col];
            if (fabs(factor) <= PIXIE_EPS) {
                continue;
            }
            for (j = 0; j <= total_columns; ++j) {
                sm->table[i][j] -= factor * sm->table[pivot_row][j];
            }
        }
    }
}

static bool pixie_simplex_is_safe(const PixieSimplex *sm) {
    return sm->table[sm->e_cnt][0] >= -PIXIE_FEAS_TOL;
}

static PixieStatus pixie_solve_lp_relaxation(const PixieModel *m, const double *node_lb, const double *node_ub, double *x_out, double *obj_out_min) {
    int n = m->n_vars;
    int i;
    int k;
    int ny = 0;
    double *lb_used = NULL;
    double *ub_used = NULL;
    PixieVarMap *map = NULL;
    PixieConVec tcons;
    PixieConstraint tc;
    double *obj_y = NULL;
    double *obj_max = NULL;
    PixieSimplex sm;
    PixieSimplexStatus ss;
    double *y = NULL;
    PixieStatus ret = PIXIE_STATUS_ERROR;

    lb_used = (double *)pixie_xmalloc((size_t)n * sizeof(double));
    ub_used = (double *)pixie_xmalloc((size_t)n * sizeof(double));
    map = (PixieVarMap *)pixie_xmalloc((size_t)n * sizeof(PixieVarMap));
    pixie_convec_init(&tcons);
    memset(&sm, 0, sizeof(sm));

    for (i = 0; i < n; ++i) {
        double lb = m->vars[i].lb;
        double ub = m->vars[i].ub;
        if (node_lb != NULL) {
            lb = pixie_max(lb, node_lb[i]);
        }
        if (node_ub != NULL) {
            ub = pixie_min(ub, node_ub[i]);
        }
        if (m->vars[i].type == PIXIE_VAR_BIN) {
            lb = pixie_max(lb, 0.0);
            ub = pixie_min(ub, 1.0);
        }
        if (m->vars[i].type != PIXIE_VAR_CONT) {
            if (pixie_is_finite(lb)) {
                lb = ceil(lb - PIXIE_INT_TOL);
            }
            if (pixie_is_finite(ub)) {
                ub = floor(ub + PIXIE_INT_TOL);
            }
        }
        if (pixie_is_finite(lb) && pixie_is_finite(ub) && lb > ub + PIXIE_FEAS_TOL) {
            ret = PIXIE_STATUS_INFEASIBLE;
            goto cleanup;
        }
        lb_used[i] = lb;
        ub_used[i] = ub;
        if (pixie_is_finite(lb)) {
            map[i].kind = PIXIE_MAP_SHIFT_LB;
            map[i].y0 = ny++;
            map[i].y1 = -1;
            map[i].constant = lb;
        } else if (pixie_is_finite(ub)) {
            map[i].kind = PIXIE_MAP_SHIFT_UB;
            map[i].y0 = ny++;
            map[i].y1 = -1;
            map[i].constant = ub;
        } else {
            map[i].kind = PIXIE_MAP_FREE;
            map[i].y0 = ny++;
            map[i].y1 = ny++;
            map[i].constant = 0.0;
        }
    }

    obj_y = (double *)pixie_xcalloc((size_t)ny, sizeof(double));
    obj_max = (double *)pixie_xcalloc((size_t)ny, sizeof(double));

    for (i = 0; i < n; ++i) {
        double cmin = (m->obj_sense == +1) ? m->vars[i].obj : -m->vars[i].obj;
        if (map[i].kind == PIXIE_MAP_SHIFT_LB) {
            obj_y[map[i].y0] += cmin;
        } else if (map[i].kind == PIXIE_MAP_SHIFT_UB) {
            obj_y[map[i].y0] += -cmin;
        } else {
            obj_y[map[i].y0] += cmin;
            obj_y[map[i].y1] += -cmin;
        }
    }
    for (i = 0; i < ny; ++i) {
        obj_max[i] = -obj_y[i];
    }

    for (i = 0; i < m->n_cons; ++i) {
        tc.cmp = m->cons[i].cmp;
        tc.rhs = m->cons[i].rhs;
        pixie_sparse_init(&tc.a);
        for (k = 0; k < m->cons[i].a.len; ++k) {
            int v = m->cons[i].a.idx[k];
            double a = m->cons[i].a.val[k];
            if (map[v].kind == PIXIE_MAP_SHIFT_LB) {
                pixie_sparse_add(&tc.a, map[v].y0, a);
                tc.rhs -= a * map[v].constant;
            } else if (map[v].kind == PIXIE_MAP_SHIFT_UB) {
                pixie_sparse_add(&tc.a, map[v].y0, -a);
                tc.rhs -= a * map[v].constant;
            } else {
                pixie_sparse_add(&tc.a, map[v].y0, a);
                pixie_sparse_add(&tc.a, map[v].y1, -a);
            }
        }
        pixie_sparse_prune(&tc.a, PIXIE_EPS);
        pixie_convec_push(&tcons, &tc);
    }

    for (i = 0; i < n; ++i) {
        if (map[i].kind == PIXIE_MAP_SHIFT_LB && pixie_is_finite(ub_used[i])) {
            double rhs = ub_used[i] - lb_used[i];
            if (rhs < -PIXIE_FEAS_TOL) {
                ret = PIXIE_STATUS_INFEASIBLE;
                goto cleanup;
            }
            tc.cmp = PIXIE_CMP_LE;
            tc.rhs = rhs;
            pixie_sparse_init(&tc.a);
            pixie_sparse_add(&tc.a, map[i].y0, 1.0);
            pixie_convec_push(&tcons, &tc);
        }
    }

    if (!pixie_make_simplex_table(tcons.data, tcons.len, ny, obj_max, &sm)) {
        ret = PIXIE_STATUS_ERROR;
        goto cleanup;
    }

    if (sm.a_cnt > 0) {
        ss = pixie_simplex_solve(&sm, ny + sm.s_cnt + sm.a_cnt);
        if (ss != PIXIE_SIMPLEX_OPTIMAL) {
            ret = PIXIE_STATUS_INFEASIBLE;
            goto cleanup;
        }
        if (!pixie_simplex_is_safe(&sm)) {
            ret = PIXIE_STATUS_INFEASIBLE;
            goto cleanup;
        }
        pixie_change_z_param(&sm, ny, obj_max);
        ss = pixie_simplex_solve(&sm, ny + sm.s_cnt);
        if (ss == PIXIE_SIMPLEX_UNBOUNDED) {
            ret = PIXIE_STATUS_UNBOUNDED;
            goto cleanup;
        }
        if (ss != PIXIE_SIMPLEX_OPTIMAL) {
            ret = PIXIE_STATUS_UNKNOWN;
            goto cleanup;
        }
    } else {
        ss = pixie_simplex_solve(&sm, ny + sm.s_cnt);
        if (ss == PIXIE_SIMPLEX_UNBOUNDED) {
            ret = PIXIE_STATUS_UNBOUNDED;
            goto cleanup;
        }
        if (ss != PIXIE_SIMPLEX_OPTIMAL) {
            ret = PIXIE_STATUS_UNKNOWN;
            goto cleanup;
        }
    }

    y = (double *)pixie_xcalloc((size_t)ny, sizeof(double));
    for (i = 0; i < sm.e_cnt; ++i) {
        if (sm.v_idx[i] >= 0 && sm.v_idx[i] < ny) {
            y[sm.v_idx[i]] = sm.table[i][0];
        }
    }

    for (i = 0; i < n; ++i) {
        double xv;
        if (map[i].kind == PIXIE_MAP_SHIFT_LB) {
            xv = map[i].constant + y[map[i].y0];
        } else if (map[i].kind == PIXIE_MAP_SHIFT_UB) {
            xv = map[i].constant - y[map[i].y0];
        } else {
            xv = y[map[i].y0] - y[map[i].y1];
        }
        if (pixie_is_finite(lb_used[i]) && xv < lb_used[i] && xv > lb_used[i] - 1e-6) {
            xv = lb_used[i];
        }
        if (pixie_is_finite(ub_used[i]) && xv > ub_used[i] && xv < ub_used[i] + 1e-6) {
            xv = ub_used[i];
        }
        x_out[i] = xv;
    }

    *obj_out_min = 0.0;
    for (i = 0; i < n; ++i) {
        double cmin = (m->obj_sense == +1) ? m->vars[i].obj : -m->vars[i].obj;
        *obj_out_min += cmin * x_out[i];
    }
    ret = PIXIE_STATUS_OPTIMAL;

cleanup:
    free(lb_used);
    free(ub_used);
    free(map);
    pixie_convec_free(&tcons);
    free(obj_y);
    free(obj_max);
    if (sm.table != NULL) {
        pixie_simplex_free(&sm);
    }
    free(y);
    return ret;
}

static void pixie_nodestack_init(PixieNodeStack *s) {
    s->data = NULL;
    s->len = 0;
    s->cap = 0;
}

static void pixie_nodestack_free(PixieNodeStack *s) {
    int i;
    for (i = 0; i < s->len; ++i) {
        free(s->data[i].lb);
        free(s->data[i].ub);
    }
    free(s->data);
    s->data = NULL;
    s->len = 0;
    s->cap = 0;
}

static void pixie_nodestack_push_copy(PixieNodeStack *s, const double *lb, const double *ub, int n, double bound, int depth) {
    int nc;
    PixieNode *nd;
    if (s->len >= s->cap) {
        nc = (s->cap <= 0) ? 8 : s->cap;
        while (nc <= s->len) {
            if (nc > INT_MAX / 2) {
                nc = s->len + 1;
                break;
            }
            nc *= 2;
        }
        s->data = (PixieNode *)pixie_xrealloc(s->data, (size_t)nc * sizeof(PixieNode));
        s->cap = nc;
    }
    nd = &s->data[s->len];
    nd->lb = (double *)pixie_xmalloc((size_t)n * sizeof(double));
    nd->ub = (double *)pixie_xmalloc((size_t)n * sizeof(double));
    memcpy(nd->lb, lb, (size_t)n * sizeof(double));
    memcpy(nd->ub, ub, (size_t)n * sizeof(double));
    nd->bound = bound;
    nd->depth = depth;
    ++s->len;
}

static PixieNode pixie_nodestack_pop(PixieNodeStack *s) {
    PixieNode nd;
    nd.lb = NULL;
    nd.ub = NULL;
    nd.bound = 0.0;
    nd.depth = 0;
    if (s->len <= 0) {
        return nd;
    }
    --s->len;
    nd = s->data[s->len];
    return nd;
}

static double pixie_nodestack_min_bound(const PixieNodeStack *s) {
    int i;
    double b = HUGE_VAL;
    for (i = 0; i < s->len; ++i) {
        if (s->data[i].bound < b) {
            b = s->data[i].bound;
        }
    }
    return b;
}

static bool pixie_is_integral(double x) {
    return fabs(x - nearbyint(x)) <= PIXIE_INT_TOL;
}

static int pixie_choose_branch_var(const PixieModel *m, const double *x, bool use_seed, uint64_t *rng_state, int *frac_cnt_out) {
    int i;
    int pick = -1;
    int frac_cnt = 0;
    int tie_cnt = 0;
    double best_score = -1.0;
    for (i = 0; i < m->n_vars; ++i) {
        double xv;
        double f;
        double score;
        if (m->vars[i].type == PIXIE_VAR_CONT) {
            continue;
        }
        xv = x[i];
        if (pixie_is_integral(xv)) {
            continue;
        }
        ++frac_cnt;
        f = fabs(xv - floor(xv));
        if (f > 0.5) {
            f = 1.0 - f;
        }
        score = f;
        if (pick < 0 || score > best_score + 1e-15) {
            pick = i;
            best_score = score;
            tie_cnt = 1;
        } else if (fabs(score - best_score) <= 1e-15) {
            if (use_seed) {
                ++tie_cnt;
                if ((pixie_rng_next(rng_state) % (uint64_t)tie_cnt) == 0ULL) {
                    pick = i;
                }
            } else if (i < pick) {
                pick = i;
            }
        }
    }
    if (frac_cnt_out != NULL) {
        *frac_cnt_out = frac_cnt;
    }
    return pick;
}

static bool pixie_try_rounding_heuristic(const PixieModel *m,
                                         const double *node_lb,
                                         const double *node_ub,
                                         const double *lp_x,
                                         double *cand_x,
                                         double *obj_min_out) {
    int i;
    int k;
    for (i = 0; i < m->n_vars; ++i) {
        double lb = m->vars[i].lb;
        double ub = m->vars[i].ub;
        double xv = lp_x[i];
        if (node_lb != NULL) {
            lb = pixie_max(lb, node_lb[i]);
        }
        if (node_ub != NULL) {
            ub = pixie_min(ub, node_ub[i]);
        }
        if (m->vars[i].type == PIXIE_VAR_BIN) {
            lb = pixie_max(lb, 0.0);
            ub = pixie_min(ub, 1.0);
        }
        if (m->vars[i].type != PIXIE_VAR_CONT) {
            if (pixie_is_finite(lb)) {
                lb = ceil(lb - PIXIE_INT_TOL);
            }
            if (pixie_is_finite(ub)) {
                ub = floor(ub + PIXIE_INT_TOL);
            }
        }
        if (pixie_is_finite(lb) && pixie_is_finite(ub) && lb > ub + PIXIE_FEAS_TOL) {
            return false;
        }
        if (pixie_is_finite(lb) && xv < lb) {
            xv = lb;
        }
        if (pixie_is_finite(ub) && xv > ub) {
            xv = ub;
        }
        if (m->vars[i].type != PIXIE_VAR_CONT) {
            xv = nearbyint(xv);
            if (pixie_is_finite(lb) && xv < lb) {
                xv = lb;
            }
            if (pixie_is_finite(ub) && xv > ub) {
                xv = ub;
            }
            if (!pixie_is_integral(xv)) {
                return false;
            }
        }
        cand_x[i] = xv;
    }

    for (i = 0; i < m->n_cons; ++i) {
        double lhs = 0.0;
        double rhs = m->cons[i].rhs;
        double tol = 1e-6 + 1e-9 * fabs(rhs);
        for (k = 0; k < m->cons[i].a.len; ++k) {
            int v = m->cons[i].a.idx[k];
            lhs += m->cons[i].a.val[k] * cand_x[v];
        }
        if (m->cons[i].cmp == PIXIE_CMP_EQ) {
            if (fabs(lhs - rhs) > tol) {
                return false;
            }
        } else if (m->cons[i].cmp == PIXIE_CMP_LE) {
            if (lhs > rhs + tol) {
                return false;
            }
        } else {
            if (lhs < rhs - tol) {
                return false;
            }
        }
    }

    *obj_min_out = 0.0;
    for (i = 0; i < m->n_vars; ++i) {
        double cmin = (m->obj_sense == +1) ? m->vars[i].obj : -m->vars[i].obj;
        *obj_min_out += cmin * cand_x[i];
    }
    return true;
}

static double pixie_obj_out_from_min(const PixieModel *m, double obj_min) {
    return (m->obj_sense == +1) ? obj_min : -obj_min;
}

static PixieSolution pixie_solution_create(int n) {
    PixieSolution sol;
    sol.status = PIXIE_STATUS_UNKNOWN;
    sol.obj_min = NAN;
    sol.obj_out = NAN;
    sol.x = (double *)pixie_xcalloc((size_t)n, sizeof(double));
    sol.n = n;
    sol.has_primal = false;
    sol.nodes_processed = 0;
    sol.stopped_time = false;
    sol.stopped_nodes = false;
    sol.stopped_gap = false;
    return sol;
}

static void pixie_solution_free(PixieSolution *sol) {
    free(sol->x);
    sol->x = NULL;
    sol->n = 0;
}

static PixieSolution pixie_solve(const PixieModel *m, const PixieOptions *opt) {
    PixieSolution sol = pixie_solution_create(m->n_vars);
    int n = m->n_vars;
    int n_int = pixie_count_integer_vars(m);
    double *lp_x = (double *)pixie_xmalloc((size_t)n * sizeof(double));
    double *heur_x = NULL;
    double *root_lb = NULL;
    double *root_ub = NULL;
    double lp_obj = NAN;
    PixieStatus lps;
    int i;

    root_lb = (double *)pixie_xmalloc((size_t)n * sizeof(double));
    root_ub = (double *)pixie_xmalloc((size_t)n * sizeof(double));
    for (i = 0; i < n; ++i) {
        root_lb[i] = m->vars[i].lb;
        root_ub[i] = m->vars[i].ub;
        if (m->vars[i].type == PIXIE_VAR_BIN) {
            root_lb[i] = pixie_max(root_lb[i], 0.0);
            root_ub[i] = pixie_min(root_ub[i], 1.0);
        }
    }

    if (opt->pure_lp || n_int == 0) {
        lps = pixie_solve_lp_relaxation(m, root_lb, root_ub, lp_x, &lp_obj);
        sol.nodes_processed = 1;
        if (lps == PIXIE_STATUS_OPTIMAL) {
            sol.status = PIXIE_STATUS_OPTIMAL;
            sol.obj_min = lp_obj;
            sol.obj_out = pixie_obj_out_from_min(m, lp_obj);
            memcpy(sol.x, lp_x, (size_t)n * sizeof(double));
            sol.has_primal = true;
        } else if (lps == PIXIE_STATUS_INFEASIBLE) {
            sol.status = PIXIE_STATUS_INFEASIBLE;
        } else if (lps == PIXIE_STATUS_UNBOUNDED) {
            sol.status = PIXIE_STATUS_UNBOUNDED;
        } else {
            sol.status = PIXIE_STATUS_UNKNOWN;
        }
        free(lp_x);
        free(root_lb);
        free(root_ub);
        return sol;
    }

    heur_x = (double *)pixie_xmalloc((size_t)n * sizeof(double));

    {
        PixieNodeStack stack;
        double best_obj = HUGE_VAL;
        bool have_inc = false;
        bool saw_unbounded_relax = false;
        uint64_t rng_state = opt->seed_set ? opt->seed : UINT64_C(0x6a09e667f3bcc909);
        clock_t t0 = clock();

        pixie_nodestack_init(&stack);
        pixie_nodestack_push_copy(&stack, root_lb, root_ub, n, -HUGE_VAL, 0);

        while (stack.len > 0) {
            PixieNode node;
            bool time_stop = false;
            bool node_stop = false;

            if (opt->node_limit > 0 && sol.nodes_processed >= opt->node_limit) {
                node_stop = true;
            }
            if (opt->time_limit_sec > 0.0) {
                double elapsed = (double)(clock() - t0) / (double)CLOCKS_PER_SEC;
                if (elapsed >= opt->time_limit_sec) {
                    time_stop = true;
                }
            }
            if (time_stop || node_stop) {
                sol.stopped_time = time_stop;
                sol.stopped_nodes = node_stop;
                break;
            }

            node = pixie_nodestack_pop(&stack);
            if (node.lb == NULL || node.ub == NULL) {
                continue;
            }

            if (have_inc && node.bound >= best_obj - PIXIE_FEAS_TOL) {
                free(node.lb);
                free(node.ub);
                continue;
            }

            lps = pixie_solve_lp_relaxation(m, node.lb, node.ub, lp_x, &lp_obj);
            ++sol.nodes_processed;
            if (opt->verbose >= 2 && (sol.nodes_processed % 1000LL) == 0LL) {
                fprintf(stderr, "c nodes=%lld incumbent=%s\n",
                        sol.nodes_processed,
                        have_inc ? "yes" : "no");
            }

            if (lps == PIXIE_STATUS_INFEASIBLE) {
                free(node.lb);
                free(node.ub);
                continue;
            }
            if (lps == PIXIE_STATUS_UNBOUNDED) {
                saw_unbounded_relax = true;
                free(node.lb);
                free(node.ub);
                continue;
            }
            if (lps != PIXIE_STATUS_OPTIMAL) {
                free(node.lb);
                free(node.ub);
                continue;
            }

            if (have_inc && lp_obj >= best_obj - PIXIE_FEAS_TOL) {
                free(node.lb);
                free(node.ub);
                continue;
            }

            {
                int frac_cnt = 0;
                int bvar = pixie_choose_branch_var(m, lp_x, opt->seed_set, &rng_state, &frac_cnt);
                if (bvar < 0) {
                    if (!have_inc || lp_obj < best_obj - PIXIE_FEAS_TOL) {
                        best_obj = lp_obj;
                        have_inc = true;
                        memcpy(sol.x, lp_x, (size_t)n * sizeof(double));
                        sol.has_primal = true;
                        if (opt->verbose >= 1) {
                            fprintf(stderr, "c incumbent %.17g at node %lld\n", pixie_obj_out_from_min(m, best_obj), sol.nodes_processed);
                        }
                    }
                } else {
                    bool pushed = false;
                    double v = lp_x[bvar];
                    double lo = floor(v);
                    double hi = ceil(v);
                    double old_lb = node.lb[bvar];
                    double old_ub = node.ub[bvar];
                    bool can_down = false;
                    bool can_up = false;
                    double ub_down = 0.0;
                    double lb_up = 0.0;

                    if (!have_inc && frac_cnt <= 24 && node.depth <= 6) {
                        double heur_obj = NAN;
                        if (pixie_try_rounding_heuristic(m, node.lb, node.ub, lp_x, heur_x, &heur_obj)) {
                            if (!have_inc || heur_obj < best_obj - PIXIE_FEAS_TOL) {
                                best_obj = heur_obj;
                                have_inc = true;
                                memcpy(sol.x, heur_x, (size_t)n * sizeof(double));
                                sol.has_primal = true;
                                if (opt->verbose >= 1) {
                                    fprintf(stderr, "c incumbent %.17g at node %lld (rounding)\n", pixie_obj_out_from_min(m, best_obj), sol.nodes_processed);
                                }
                            }
                        }
                    }

                    if (have_inc && lp_obj >= best_obj - PIXIE_FEAS_TOL) {
                        pushed = true;
                    } else {
                        if (lo >= old_lb - PIXIE_FEAS_TOL) {
                            ub_down = pixie_min(old_ub, lo);
                            if (!(pixie_is_finite(old_lb) && pixie_is_finite(ub_down) && old_lb > ub_down + PIXIE_FEAS_TOL)) {
                                can_down = true;
                            }
                        }
                        if (hi <= old_ub + PIXIE_FEAS_TOL) {
                            lb_up = pixie_max(old_lb, hi);
                            if (!(pixie_is_finite(lb_up) && pixie_is_finite(old_ub) && lb_up > old_ub + PIXIE_FEAS_TOL)) {
                                can_up = true;
                            }
                        }

                        if (can_down) {
                            node.lb[bvar] = old_lb;
                            node.ub[bvar] = ub_down;
                            pixie_nodestack_push_copy(&stack, node.lb, node.ub, n, lp_obj, node.depth + 1);
                            pushed = true;
                        }
                        if (can_up) {
                            node.lb[bvar] = lb_up;
                            node.ub[bvar] = old_ub;
                            pixie_nodestack_push_copy(&stack, node.lb, node.ub, n, lp_obj, node.depth + 1);
                            pushed = true;
                        }
                    }

                    node.lb[bvar] = old_lb;
                    node.ub[bvar] = old_ub;

                    if (!pushed) {
                        if (!have_inc || lp_obj < best_obj - PIXIE_FEAS_TOL) {
                            best_obj = lp_obj;
                            have_inc = true;
                            memcpy(sol.x, lp_x, (size_t)n * sizeof(double));
                            sol.has_primal = true;
                            if (opt->verbose >= 1) {
                                fprintf(stderr, "c incumbent %.17g at node %lld\n", pixie_obj_out_from_min(m, best_obj), sol.nodes_processed);
                            }
                        }
                    }
                }
            }

            free(node.lb);
            free(node.ub);

            if (have_inc && opt->gap_limit >= 0.0 && stack.len > 0) {
                double global_lb = pixie_nodestack_min_bound(&stack);
                double denom = fabs(best_obj) + 1e-12;
                double gap = HUGE_VAL;
                if (pixie_is_finite(global_lb)) {
                    if (global_lb > best_obj) {
                        global_lb = best_obj;
                    }
                    gap = (best_obj - global_lb) / denom;
                }
                if (gap <= opt->gap_limit + 1e-15) {
                    sol.stopped_gap = true;
                    break;
                }
            }
        }

        if (have_inc) {
            sol.obj_min = best_obj;
            sol.obj_out = pixie_obj_out_from_min(m, best_obj);
            if (stack.len == 0 && !sol.stopped_time && !sol.stopped_nodes && !sol.stopped_gap && !saw_unbounded_relax) {
                sol.status = PIXIE_STATUS_OPTIMAL;
            } else {
                sol.status = PIXIE_STATUS_UNKNOWN;
            }
        } else {
            if (sol.stopped_time || sol.stopped_nodes || sol.stopped_gap) {
                sol.status = PIXIE_STATUS_UNKNOWN;
            } else if (stack.len == 0 && !saw_unbounded_relax) {
                sol.status = PIXIE_STATUS_INFEASIBLE;
            } else {
                sol.status = PIXIE_STATUS_UNKNOWN;
            }
        }

        pixie_nodestack_free(&stack);
    }

    free(lp_x);
    free(heur_x);
    free(root_lb);
    free(root_ub);
    return sol;
}

static void pixie_print_solution(const PixieModel *m, const PixieSolution *sol) {
    int i;
    (void)m;
    if (sol->status == PIXIE_STATUS_OPTIMAL) {
        printf("s OPTIMUM FOUND\n");
        printf("o %.17g\n", sol->obj_out);
        printf("v");
        for (i = 0; i < sol->n; ++i) {
            printf(" %.17g", sol->x[i]);
        }
        printf("\n");
        return;
    }
    if (sol->status == PIXIE_STATUS_INFEASIBLE) {
        printf("s INFEASIBLE\n");
        return;
    }
    if (sol->status == PIXIE_STATUS_UNBOUNDED) {
        printf("s UNBOUNDED\n");
        return;
    }
    printf("s UNKNOWN\n");
    if (sol->has_primal) {
        printf("o %.17g\n", sol->obj_out);
    } else {
        printf("o nan\n");
    }
}

static int pixie_selftest_one_lp(void) {
    PixieModel m;
    PixieSparse a;
    PixieOptions opt;
    PixieSolution sol;
    int x;
    int y;

    pixie_model_init(&m);
    x = pixie_model_get_var(&m, "x", true);
    y = pixie_model_get_var(&m, "y", true);
    m.obj_sense = +1;
    m.vars[x].obj = 1.0;
    m.vars[y].obj = 1.0;
    pixie_sparse_init(&a);
    pixie_sparse_add(&a, x, 1.0);
    pixie_sparse_add(&a, y, 1.0);
    pixie_model_add_constraint(&m, &a, PIXIE_CMP_GE, 1.0);
    pixie_sparse_free(&a);

    memset(&opt, 0, sizeof(opt));
    opt.gap_limit = -1.0;
    sol = pixie_solve(&m, &opt);
    if (sol.status != PIXIE_STATUS_OPTIMAL || fabs(sol.obj_out - 1.0) > 1e-6) {
        pixie_solution_free(&sol);
        pixie_model_free(&m);
        return 0;
    }
    pixie_solution_free(&sol);
    pixie_model_free(&m);
    return 1;
}

static int pixie_selftest_one_mip(void) {
    PixieModel m;
    PixieSparse a;
    PixieOptions opt;
    PixieSolution sol;
    int x;
    int y;

    pixie_model_init(&m);
    x = pixie_model_get_var(&m, "x", true);
    y = pixie_model_get_var(&m, "y", true);
    m.obj_sense = +1;
    m.vars[x].obj = 1.0;
    m.vars[y].obj = 1.0;
    pixie_var_mark_binary(&m, x);
    pixie_var_mark_binary(&m, y);
    pixie_sparse_init(&a);
    pixie_sparse_add(&a, x, 1.0);
    pixie_sparse_add(&a, y, 1.0);
    pixie_model_add_constraint(&m, &a, PIXIE_CMP_GE, 1.0);
    pixie_sparse_free(&a);

    memset(&opt, 0, sizeof(opt));
    opt.gap_limit = -1.0;
    sol = pixie_solve(&m, &opt);
    if ((sol.status != PIXIE_STATUS_OPTIMAL && sol.status != PIXIE_STATUS_UNKNOWN) || !sol.has_primal || fabs(sol.obj_out - 1.0) > 1e-6) {
        pixie_solution_free(&sol);
        pixie_model_free(&m);
        return 0;
    }
    pixie_solution_free(&sol);
    pixie_model_free(&m);
    return 1;
}

static int pixie_selftest_infeasible(void) {
    PixieModel m;
    PixieSparse a;
    PixieOptions opt;
    PixieSolution sol;
    int x;

    pixie_model_init(&m);
    x = pixie_model_get_var(&m, "x", true);
    m.obj_sense = +1;
    m.vars[x].obj = 1.0;
    pixie_sparse_init(&a);
    pixie_sparse_add(&a, x, 1.0);
    pixie_model_add_constraint(&m, &a, PIXIE_CMP_LE, -1.0);
    pixie_sparse_free(&a);

    memset(&opt, 0, sizeof(opt));
    opt.gap_limit = -1.0;
    sol = pixie_solve(&m, &opt);
    if (sol.status != PIXIE_STATUS_INFEASIBLE) {
        pixie_solution_free(&sol);
        pixie_model_free(&m);
        return 0;
    }
    pixie_solution_free(&sol);
    pixie_model_free(&m);
    return 1;
}

static int pixie_selftest_unbounded(void) {
    PixieModel m;
    PixieOptions opt;
    PixieSolution sol;
    int x;

    pixie_model_init(&m);
    x = pixie_model_get_var(&m, "x", true);
    m.obj_sense = +1;
    m.vars[x].obj = -1.0;

    memset(&opt, 0, sizeof(opt));
    opt.pure_lp = true;
    opt.gap_limit = -1.0;
    sol = pixie_solve(&m, &opt);
    if (sol.status != PIXIE_STATUS_UNBOUNDED) {
        pixie_solution_free(&sol);
        pixie_model_free(&m);
        return 0;
    }
    pixie_solution_free(&sol);
    pixie_model_free(&m);
    return 1;
}

static int pixie_run_selftest(void) {
    if (!pixie_selftest_one_lp()) {
        fprintf(stderr, "selftest failed: lp optimum\n");
        return 1;
    }
    if (!pixie_selftest_one_mip()) {
        fprintf(stderr, "selftest failed: mip optimum\n");
        return 1;
    }
    if (!pixie_selftest_infeasible()) {
        fprintf(stderr, "selftest failed: infeasible\n");
        return 1;
    }
    if (!pixie_selftest_unbounded()) {
        fprintf(stderr, "selftest failed: unbounded\n");
        return 1;
    }
    fprintf(stderr, "selftest: OK\n");
    return 0;
}

static void pixie_print_usage(const char *prog) {
    fprintf(stderr,
            "usage:\n"
            "  %s <file>\n"
            "  %s --lp <file>\n"
            "  %s --mps <file>\n"
            "options:\n"
            "  --time <sec>     time limit\n"
            "  --node <n>       node limit\n"
            "  --gap <g>        relative MIP gap limit\n"
            "  --seed <s>       RNG seed (enables randomized branching tie-breaks)\n"
            "  --verbose <k>    verbosity 0..3\n"
            "  --purelp         solve LP relaxation only\n"
            "  --selftest       run embedded tests\n",
            prog, prog, prog);
}

int main(int argc, char **argv) {
    PixieOptions opt;
    PixieModel model;
    PixieSolution sol;
    char err[512];
    int i;

    memset(&opt, 0, sizeof(opt));
    opt.format = PIXIE_FORMAT_AUTO;
    opt.file = NULL;
    opt.time_limit_sec = 0.0;
    opt.node_limit = 0;
    opt.gap_limit = -1.0;
    opt.seed = 0;
    opt.seed_set = false;
    opt.verbose = 0;
    opt.pure_lp = false;
    opt.selftest = false;

    for (i = 1; i < argc; ++i) {
        const char *a = argv[i];
        if (strcmp(a, "--selftest") == 0) {
            opt.selftest = true;
        } else if (strcmp(a, "--lp") == 0) {
            if (i + 1 >= argc) {
                pixie_print_usage(argv[0]);
                return EXIT_FAILURE;
            }
            opt.format = PIXIE_FORMAT_LP;
            opt.file = argv[++i];
        } else if (strcmp(a, "--mps") == 0) {
            if (i + 1 >= argc) {
                pixie_print_usage(argv[0]);
                return EXIT_FAILURE;
            }
            opt.format = PIXIE_FORMAT_MPS;
            opt.file = argv[++i];
        } else if (strcmp(a, "--time") == 0) {
            if (i + 1 >= argc || !pixie_parse_double_token(argv[i + 1], &opt.time_limit_sec)) {
                fprintf(stderr, "invalid --time value\n");
                return EXIT_FAILURE;
            }
            ++i;
        } else if (strcmp(a, "--node") == 0) {
            if (i + 1 >= argc || !pixie_parse_ll_token(argv[i + 1], &opt.node_limit)) {
                fprintf(stderr, "invalid --node value\n");
                return EXIT_FAILURE;
            }
            ++i;
        } else if (strcmp(a, "--gap") == 0) {
            if (i + 1 >= argc || !pixie_parse_double_token(argv[i + 1], &opt.gap_limit)) {
                fprintf(stderr, "invalid --gap value\n");
                return EXIT_FAILURE;
            }
            ++i;
        } else if (strcmp(a, "--seed") == 0) {
            if (i + 1 >= argc || !pixie_parse_u64_token(argv[i + 1], &opt.seed)) {
                fprintf(stderr, "invalid --seed value\n");
                return EXIT_FAILURE;
            }
            opt.seed_set = true;
            ++i;
        } else if (strcmp(a, "--verbose") == 0) {
            long long vv = 0;
            if (i + 1 >= argc || !pixie_parse_ll_token(argv[i + 1], &vv)) {
                fprintf(stderr, "invalid --verbose value\n");
                return EXIT_FAILURE;
            }
            if (vv < 0) {
                vv = 0;
            }
            if (vv > 3) {
                vv = 3;
            }
            opt.verbose = (int)vv;
            ++i;
        } else if (strcmp(a, "--purelp") == 0) {
            opt.pure_lp = true;
        } else if (a[0] == '-') {
            fprintf(stderr, "unknown option: %s\n", a);
            pixie_print_usage(argv[0]);
            return EXIT_FAILURE;
        } else {
            if (opt.file != NULL) {
                fprintf(stderr, "multiple input files provided\n");
                pixie_print_usage(argv[0]);
                return EXIT_FAILURE;
            }
            opt.file = a;
        }
    }

    if (opt.selftest) {
        return pixie_run_selftest();
    }

    if (opt.file == NULL) {
        pixie_print_usage(argv[0]);
        return EXIT_FAILURE;
    }

    pixie_model_init(&model);
    err[0] = '\0';
    if (!pixie_read_model(opt.file, opt.format, &model, err, sizeof(err))) {
        fprintf(stderr, "error: %s\n", (err[0] != '\0') ? err : "failed to read model");
        pixie_model_free(&model);
        return EXIT_FAILURE;
    }

    for (i = 0; i < model.n_vars; ++i) {
        if (model.vars[i].type == PIXIE_VAR_BIN) {
            if (model.vars[i].lb < 0.0) {
                model.vars[i].lb = 0.0;
            }
            if (model.vars[i].ub > 1.0) {
                model.vars[i].ub = 1.0;
            }
        }
    }

    sol = pixie_solve(&model, &opt);
    pixie_print_solution(&model, &sol);
    pixie_solution_free(&sol);
    pixie_model_free(&model);
    return EXIT_SUCCESS;
}