"""Small primal-dual active-set solver for strictly convex inequality QPs.

Solves min .5 v'Hv + c'v subject to A v >= b using the KKT sign convention
H v + c - A' mu = 0, mu >= 0.
"""
import numpy as np


def solve_active_set(H, c, A, b, active_init=None, tol=1e-10, max_iter=100):
    H = np.asarray(H, float); c = np.asarray(c, float)
    A = np.asarray(A, float); b = np.asarray(b, float)
    n, m = H.shape[0], A.shape[0]
    active = set([] if active_init is None else active_init)
    history = []
    v = np.linalg.solve(H, -c)
    mu = np.zeros(m)
    for it in range(1, max_iter + 1):
        inds = sorted(active)
        if inds:
            K = np.block([[H, -A[inds].T], [A[inds], np.zeros((len(inds), len(inds)))]])
            rhs = np.r_[-c, b[inds]]
            try:
                sol = np.linalg.solve(K, rhs)
            except np.linalg.LinAlgError:
                # A dependent active set is not useful; drop its last member.
                active.remove(inds[-1]); history.append(("remove_singular", inds[-1])); continue
            v = sol[:n]; mu[:] = 0.; mu[inds] = sol[n:]
        else:
            v = np.linalg.solve(H, -c); mu[:] = 0.
        residual = A @ v - b
        bad_inactive = [(residual[j], j) for j in range(m) if j not in active and residual[j] < -tol]
        bad_active = [(mu[j], j) for j in active if mu[j] < -tol]
        history.append(("solve", tuple(sorted(active)), float(residual.min()) if m else 0.,
                        float(mu.min()) if m else 0.))
        if not bad_inactive and not bad_active:
            return v, mu, sorted(active), it, history
        if bad_active:
            # Most negative multiplier must leave the active set.
            _, j = min(bad_active)
            active.remove(j); history.append(("remove", j))
        else:
            # Most violated inactive constraint enters.
            _, j = min(bad_inactive)
            active.add(j); history.append(("add", j))
    raise RuntimeError("PDAS did not converge")


def kkt_residual(H, c, A, b, v, mu):
    stationarity = H @ v + c - A.T @ mu
    primal = np.maximum(0., b - A @ v)
    dual = np.maximum(0., -mu)
    complementarity = np.abs(mu * (A @ v - b))
    return max(np.max(np.abs(stationarity)), np.max(primal), np.max(dual), np.max(complementarity))