Linear Programming

import numpy as np
from scipy.optimize import linprog

# ── CAMPAIGN SPENDING EXAMPLE ── minimise total spending ─────────────
# Negate for maximise convention; linprog minimises by default
# Primal: minimise x1+x2+x3+x4
# Constraints (negated to <=):  -(-2x1+8x2+10x4) <= -50000  etc.

c = [1, 1, 1, 1]     # minimise x1+x2+x3+x4

# Constraints Ax <= b  (negate >= constraints)
A = [[-(-2), -8,  0, -10],    # Urban: -(votes) <= -50000
     [-5,    -2,  0,   0],    # Suburban
     [-3,     5, -10,  2]]    # Rural
b = [-50000, -100000, -25000]

result = linprog(c, A_ub=A, b_ub=b, bounds=[(0,None)]*4,
                 method='highs')
print(f"Optimal spend: ${result.fun:,.0f}")
print(f"x = {result.x.round(1)}")
# Optimal spend: ~$27,928,  x ≈ [18468, 3829, 0, 5631]


# ── LP DUALITY DEMONSTRATION ── weak duality certificate ─────────────
def verify_weak_duality(c, A, b, x_primal, y_dual):
    """
    Weak duality: c^T x <= b^T y for any feasible primal x, dual y.
    At optimality: c^T x* = b^T y* (strong duality).
    """
    primal_val = np.dot(c, x_primal)
    dual_val   = np.dot(b, y_dual)
    print(f"Primal value c^T x  = {primal_val:.4f}")
    print(f"Dual   value b^T y  = {dual_val:.4f}")
    print(f"Weak duality holds: {primal_val <= dual_val + 1e-9}")
    print(f"Strong duality (at optimum): {abs(primal_val-dual_val) < 1e-6}")


# ── MAX-FLOW AS LP ── demonstrating LP subsumes flow ─────────────────
def max_flow_lp(n, edges, source, sink, capacities):
    """
    Solve max flow as an LP (educational; scipy linprog is not sparse-aware).
    Maximise flow out of source subject to:
      - Capacity: f(e) <= c(e) for each edge
      - Conservation: inflow = outflow at each internal vertex
      - Non-negativity: f(e) >= 0
    """
    m = len(edges)
    # Variables: f[0..m-1] = flow on each edge
    # Objective: -sum of flows out of source (negate to minimise)
    c_obj = [0.0] * m
    for i, (u,v) in enumerate(edges):
        if u == source:
            c_obj[i] = -1.0   # maximise source outflow

    # Capacity constraints: f[i] <= cap[i]
    A_ub = np.eye(m)
    b_ub = capacities

    # Conservation constraints (equality): for each internal node
    A_eq, b_eq = [], []
    internal = [v for v in range(n) if v != source and v != sink]
    for node in internal:
        row = [0.0] * m
        for i, (u,v) in enumerate(edges):
            if v == node: row[i] =  1.0   # inflow
            if u == node: row[i] = -1.0   # outflow
        A_eq.append(row)
        b_eq.append(0.0)

    result = linprog(c_obj, A_ub=A_ub, b_ub=b_ub,
                     A_eq=A_eq if A_eq else None,
                     b_eq=b_eq if b_eq else None,
                     bounds=[(0,None)]*m, method='highs')
    return -result.fun if result.success else None


# ── SHORTEST PATHS AS LP ── difference constraints ────────────────────
def shortest_paths_lp(n, edges, source):
    """
    Maximise sum(d) subject to d[v] - d[u] <= w(u,v) and d[source]=0.
    Equivalent to finding shortest paths from source.
    """
    m = len(edges)
    # Variables: d[0..n-1]
    c_obj = [-1.0] * n     # maximise sum d[v]

    A_ub, b_ub = [], []
    for u, v, w in edges:
        row = [0.0]*n
        row[v] =  1.0
        row[u] = -1.0
        A_ub.append(row)
        b_ub.append(w)

    # d[source] = 0
    A_eq = [[0.0]*n]; A_eq[0][source] = 1.0
    b_eq = [0.0]

    result = linprog(c_obj, A_ub=A_ub, b_ub=b_ub,
                     A_eq=A_eq, b_eq=b_eq,
                     bounds=[(None,None)]*n, method='highs')
    return result.x if result.success else None

import numpy as np

def simplex(c, A, b):
    """
    Simplex algorithm in slack form. Maximise c^T x s.t. Ax <= b, x >= 0.
    Educational implementation — not numerically robust; use scipy for production.

    Returns (optimal_value, optimal_x) or raises if infeasible/unbounded.

    Algorithm:
    1. Introduce slack variables: B = A|I (n x (d+n)), extended cost c_ext.
    2. Pivot: find entering variable (most positive objective coefficient),
       find leaving variable (min-ratio test to stay feasible),
       update tableau.
    3. Stop when no positive objective coefficient remains.
    """
    d, n = len(c), len(b)
    # Build initial tableau: [A | I | b] with objective row [-c | 0 | 0]
    # Rows 0..n-1: constraints. Row n: objective.
    tab = np.zeros((n+1, d+n+1))
    tab[:n, :d] = A                        # constraint coefficients
    tab[:n, d:d+n] = np.eye(n)            # slack variables
    tab[:n, -1] = b                        # RHS
    tab[n, :d] = -np.array(c, dtype=float) # objective (negated for min form)

    basis = list(range(d, d+n))            # initial basis: slack variables

    for _ in range(1000):                  # max iterations guard
        # Find entering variable: most negative in objective row
        pivot_col = np.argmin(tab[n, :-1])
        if tab[n, pivot_col] >= -1e-9:
            break                          # optimal: all coefficients >= 0

        # Min-ratio test: find leaving variable
        ratios = np.full(n, np.inf)
        for i in range(n):
            if tab[i, pivot_col] > 1e-9:
                ratios[i] = tab[i, -1] / tab[i, pivot_col]
        if np.all(np.isinf(ratios)):
            raise ValueError("LP is unbounded")
        pivot_row = np.argmin(ratios)

        # Pivot: divide pivot row by pivot element, then eliminate column
        tab[pivot_row] /= tab[pivot_row, pivot_col]
        for i in range(n+1):
            if i != pivot_row:
                tab[i] -= tab[i, pivot_col] * tab[pivot_row]

        basis[pivot_row] = pivot_col

    # Extract solution
    x = np.zeros(d + n)
    for i, b_idx in enumerate(basis):
        x[b_idx] = tab[i, -1]
    return tab[n, -1], x[:d]   # objective value, original variables


# Verify on the lecture example
c_ex = [3, 1, 2]
A_ex = [[1,1,3],[2,2,5],[4,1,2]]
b_ex = [30, 24, 36]
val, x_opt = simplex(c_ex, A_ex, b_ex)
print(f"Optimal value: {val:.4f}")   # 28.0
print(f"Optimal x: {x_opt}")         # [8, 4, 0]

import numpy as np
from scipy.optimize import linprog, minimize

# ── MINIMUM COST FLOW AS LP ── ────────────────────────────────────────
def min_cost_flow_lp(n, edges, source, sink, demand, costs, capacities):
    """
    Minimum cost max-flow as LP.
    edges: list of (u, v). costs[i], capacities[i] for each edge.
    Minimise sum cost[i]*f[i] subject to:
      - Capacity: 0 <= f[i] <= capacities[i]
      - Conservation: inflow - outflow = demand at each node
    """
    m = len(edges)
    # Objective: minimise cost^T f
    c_obj = costs

    # Conservation (equality): B f = demand
    B = np.zeros((n, m))
    for i, (u, v) in enumerate(edges):
        B[u][i] = -1   # outflow
        B[v][i] =  1   # inflow

    bounds = [(0, cap) for cap in capacities]
    result = linprog(c_obj, A_eq=B, b_eq=demand,
                     bounds=bounds, method='highs')
    return result.x if result.success else None


# ── MARKOWITZ MEAN-VARIANCE (LP-RELAXED) ── ───────────────────────────
def markowitz_max_return(mu, Sigma, risk_limit, long_only=True):
    """
    Maximise expected return mu^T w subject to:
      - Portfolio risk w^T Sigma w <= risk_limit^2  (quadratic constraint!)
      - Sum weights = 1
      - w >= 0 (long only)

    Note: this is a QP (quadratic programming), not a pure LP.
    LP approximation: linearise variance using diagonal only.
    True QP: use scipy minimize with method='SLSQP'.
    """
    n = len(mu)

    def neg_return(w): return -np.dot(mu, w)
    def grad_return(w): return -mu

    constraints = [
        {'type': 'eq',  'fun': lambda w: np.sum(w) - 1},            # weights sum to 1
        {'type': 'ineq','fun': lambda w: risk_limit**2 - w@Sigma@w} # risk constraint
    ]
    bounds = [(0,1)]*n if long_only else [(None,None)]*n
    w0 = np.ones(n)/n

    result = minimize(neg_return, w0, jac=grad_return,
                      method='SLSQP', bounds=bounds,
                      constraints=constraints)
    return result.x if result.success else None


# Example: 4 NSE instruments
mu = np.array([0.12, 0.18, 0.10, 0.15])    # expected annual returns
Sigma = np.array([                           # covariance matrix
    [0.04, 0.02, 0.01, 0.02],
    [0.02, 0.09, 0.02, 0.03],
    [0.01, 0.02, 0.03, 0.01],
    [0.02, 0.03, 0.01, 0.06],
])
w_opt = markowitz_max_return(mu, Sigma, risk_limit=0.15)
print("Optimal portfolio weights:", w_opt.round(3))
print("Expected return:", np.dot(mu, w_opt).round(4))