Core estimation

Estimator

MicroPURC estimator: the fixed-point map \(\Phi\) and its iteration.

class micropurc.core.estimator.MicroPURCEstimator(network, forward_solver, estimation_config=None, active_set_config=None, projection_config=None)[source]

Bases: object

Estimator for the microPURC model with quadratic perturbation.

Each trip contributes a projected least-squares block, and the estimate is a fixed point of

\[\Phi_N(\beta) = -H_N(\beta)^{-1} g_N(\beta), \qquad H_N = \sum_i Z_i^{\top} \hat P_i Z_i, \qquad g_N = \sum_i Z_i^{\top} \hat P_i (\mathrm{scale} \odot y_i),\]

where the sums run over the trips whose active set survives grounding, \(Z_i\) is the design restricted to trip \(i\)’s active links, and \(\hat P_i\) projects onto the null space of that trip’s grounded active incidence. The dependence of \(\Phi_N\) on \(\beta\) runs through the active sets, which come from the forward solve, so \(\Phi_N\) is piecewise constant.

fit() drives \(\beta = \Phi_N(\beta)\) by damped Krasnoselskii-Mann iteration: the step is \(\beta + \lambda (\Phi_N(\beta) - \beta)\), with \(\lambda\) backtracked by a factor ls_tau until the merit \(Q(\beta) = \tfrac12 \lVert \Phi_N(\beta) - \beta \rVert^2\) drops to or below the largest merit over the last nm_window iterations. Once the backtracking budget is spent, the smallest trial step is taken as is.

Parameters:
__init__(network, forward_solver, estimation_config=None, active_set_config=None, projection_config=None)[source]

Bind the network and forward solver and precompute native arrays.

Parameters:

  • network (Network) – The network the model is estimated on.

  • forward_solver (ForwardSolver) – Solver for the inner perturbed-utility flow problem.

  • estimation_config (EstimationConfig | None) – Iteration and tolerance settings; defaults used if omitted.

  • active_set_config (ActiveSetConfig | None) – Active-set detection settings; defaults if omitted.

  • projection_config (ProjectionConfig | None) – Range-space projection settings; defaults if omitted.

Return type:

None

phi(beta, y, b=None, *, origins=None, dests=None, lazy_design=None)[source]

Evaluate the fixed-point map \(\Phi_N(\beta)\).

Parameters:

  • beta (ndarray) – Parameter vector to evaluate at, shape (K,).

  • y – Route indicators, dense (N, L) or sparse CSR (N, L).

  • b – Node imbalances, shape (N, V). Needed when origins and dests are absent.

  • origins (ndarray | None) – Origin node indices, shape (N,).

  • dests (ndarray | None) – Destination node indices, shape (N,).

  • lazy_design (LazyDesign | None) – A lazy design evaluated on the fly by the native engine.

Returns:

\(\Phi_N(\beta) = -H_N(\beta)^{-1} g_N(\beta)\), shape (K,).

Return type:

ndarray

fit(y, b=None, beta_init=None, *, origins=None, dests=None, z_variants=None, trip_variant=None, lazy_design=None)[source]

Solve \(\beta = \Phi_N(\beta)\) from beta_init.

The iteration runs until the fixed-point residual meets fp_tol, optionally scaled by \(\sqrt{K}\), or until max_iterations is spent. Convergence additionally requires the condition number of the accumulated Hessian to stay within max_cond_H, which rules out a residual that goes small only because \(\Phi_N\) has flattened at a degenerate iterate. The final estimate is accumulated once more to obtain the scores and the covariance.

Parameters:

  • y – Route indicators, dense (N, L) or sparse CSR (N, L).

  • b – Node imbalances, shape (N, V). Needed when origins and dests are absent.

  • beta_init (Sequence[float] | None) – Initial parameter guess, shape (K,); defaults to ones.

  • origins (ndarray | None) – Origin node indices, shape (N,).

  • dests (ndarray | None) – Destination node indices, shape (N,).

  • z_variants (ndarray | None) – Per-variant Z matrices, shape (n_variants, L, K).

  • trip_variant (ndarray | None) – Per-trip variant index, shape (N,). Required with z_variants.

  • lazy_design (LazyDesign | None) – A lazy design whose per-variant regressors are evaluated on the fly, keeping the (N, L, K) design out of memory. Mutually exclusive with z_variants.

Returns:

A dict holding

  • beta_hat: the final iterate, shape (K,);

  • se_hat, cov: the i.i.d. sandwich standard errors and covariance, None when the final accumulation yields no scores or a singular Hessian;

  • u_scores: per-trip scores, shape (N, K), aligned with the input trips, a dropped trip carrying a zero row;

  • kept_mask: length-N boolean, False at a dropped trip;

  • hessian: the accumulated Hessian, shape (K, K);

  • diagnostics: the per-run summary with the convergence verdict, the merit history, and the formatted convergence and active-set tables.

Clustered inference follows from passing hessian and u_scores to cluster_robust_covariance() with length-N cluster labels: the zero rows keep dropped trips out of the meat, so the caller needs no drop bookkeeping. u_scores and kept_mask are None when the final accumulation returned no scores; if the survivor mask cannot be matched to those scores, u_scores keeps only the surviving rows and kept_mask stays None.

Raises:

ValueError – If z_variants is given without trip_variant, or if both lazy_design and z_variants are given.

Return type:

dict[str, Any]

Projection

One trip’s contribution to the estimating equations, via a Schur complement.

Every contribution is a contraction with the range-space projector

\[\hat P = I - \tilde A^{\top} \big(\tilde A \tilde A^{\top}\big)^{-1} \tilde A ,\]

where \(\tilde A\) is the trip’s active incidence submatrix after one node row per connected component has been grounded. The C++ extension applies \(\hat P\) through a Cholesky factor of the \(m \times m\) Gram matrix \(\tilde A \tilde A^{\top}\), so the \(e \times e\) projector itself is never formed. This module adapts numpy arrays to that call and names the three outputs it returns.

class micropurc.core.projection.TripContribution(H_i, g_i, u_i)[source]

Bases: object

Per-trip contribution to the estimator.

Variables:

  • H_i (numpy.ndarray) – Hessian block \(Z_{\mathrm{act}}^{\top} \hat P Z_{\mathrm{act}}\), shape (K, K).

  • g_i (numpy.ndarray) – Gradient block \(Z_{\mathrm{act}}^{\top} \hat P (\mathrm{scale} \odot y)\), shape (K,).

  • u_i (numpy.ndarray | None) – Score \(Z_{\mathrm{act}}^{\top} \hat P \big(\mathrm{scale} \odot (y - \hat x)\big)\), shape (K,); None when the flows were withheld.

Parameters:
H_i: ndarray
g_i: ndarray
u_i: ndarray | None
micropurc.core.projection.trip_contribution_schur(Z_act, A_local, y_act, scale_act, x_hat_act=None)[source]

Compute one trip’s (H_i, g_i, u_i) via the native Schur solve.

Parameters:

  • Z_act (ndarray) – Link attributes on active edges, shape (e, K).

  • A_local (ndarray) – Active incidence submatrix, shape (m, e), dense, with one node row per connected component already dropped so that \(\tilde A \tilde A^{\top}\) is invertible.

  • y_act (ndarray) – Route indicator on active edges, shape (e,).

  • scale_act (ndarray) – Perturbation scale on active edges, shape (e,).

  • x_hat_act (ndarray | None) – Optimal flows on active edges, shape (e,). Supplying it requests the score u_i; None leaves u_i unset and skips the extra right-hand side in the Schur solve.

Returns:

The trip’s H_i, g_i, and u_i.

Return type:

TripContribution

Covariance

Sandwich covariance estimation and cluster-robust standard errors.

micropurc.core.covariance.sandwich_covariance(H, u_vectors)[source]

Compute the sandwich variance and its standard errors.

\[\mathrm{cov} = H^{-1} S H^{-1}, \qquad S = \sum_i u_i u_i^{\top},\]

where each per-trip score \(u_i = Z_{\mathrm{act}}^{\top} \hat P \big(\mathrm{scale} \odot (y_i - \hat x_i)\big)\) is a \(K\)-vector, so the meat accumulates as a sum of rank-one terms.

Parameters:

  • H (ndarray) – Accumulated Hessian \(\sum_i H_i\), shape (K, K).

  • u_vectors (ndarray) – Stacked score vectors \(u_i\), shape (N, K).

Returns:

Covariance matrix, shape (K, K). se: Standard errors, shape (K,).

Return type:

cov

micropurc.core.covariance.cluster_robust_covariance(H, u_vectors, cluster_ids)[source]

Compute the cluster-robust sandwich variance and its standard errors.

Scores are summed within each cluster before the meat is formed, so within-cluster dependence is carried through to the variance:

\[S_{\mathrm{cl}} = \frac{J}{J-1} \sum_{j=1}^{J} \Big(\sum_{i \in j} u_i\Big) \Big(\sum_{i \in j} u_i\Big)^{\top}, \qquad \mathrm{cov} = H^{-1} S_{\mathrm{cl}} H^{-1},\]

over \(J\) clusters. The CR1 finite-cluster scaling \(J/(J-1)\) is included here, so a caller that applies its own would double-count it.

Both H and the meat are accumulated over observations, so the returned covariance is already scaled as a sampling variance and needs no further division by the sample size.

Wald intervals built from these standard errors take their critical value from \(t_{J-1}\).

Parameters:

  • H (ndarray) – Accumulated Hessian, shape (K, K).

  • u_vectors (ndarray) – Stacked u_i vectors, shape (N, K).

  • cluster_ids (ndarray) – Cluster (e.g., user) ID per observation, shape (N,).

Returns:

Cluster-robust covariance, shape (K, K). se: Cluster-robust standard errors, shape (K,).

Return type:

cov

Raises:

ValueError – If cluster_ids and u_vectors disagree in length, or if there are fewer than two clusters (the CR1 correction is undefined).

Active set

Active-set determination for a single trip.

A link carries flow when its dual slack is negative: with link costs \(c = Z\beta\) and the node potentials \(\lambda\) returned by the forward solve, link \(j\) is active when

\[s_j = c_j + (A^{\top}\lambda)_j < 0 .\]

The C++ extension performs the slack test, the acyclicity and destination-reachability checks on the active subnetwork, and the component-aware grounding that makes \(\tilde A \tilde A^{\top}\) invertible. This module adapts numpy and Network objects to those calls and keeps the Python result types.

class micropurc.core.activeset.ActiveSet(indices, slack, margins, kept_node_rows, num_components, valid, drop_reason='')[source]

Bases: object

Result of active-set determination for a single trip.

Variables:

  • indices (numpy.ndarray) – Active (flow-carrying) link ids, those with \(s_j < 0\).

  • slack (numpy.ndarray) – Dual slack \(s_j = c_j + (A^{\top}\lambda)_j\) for every link, shape (L,).

  • margins (numpy.ndarray | None) – Per-link margin \(|s_j| / \lVert Z_j \rVert\), shape (L,). Slack is affine in \(\beta\) through \(c = Z\beta\), so with the node potentials held fixed this is the length of the smallest parameter perturbation that flips link \(j\) between active and inactive. Links with \(\lVert Z_j \rVert = 0\) carry inf, so the trip’s margin is the minimum of the vector. None unless z_rownorm is supplied to determine_active_set().

  • kept_node_rows (numpy.ndarray) – Node rows of \(A\) retained for the range-space projection, one row per connected component having been dropped.

  • num_components (int) – Number of connected components in the active subnetwork.

  • valid (bool) – Whether the trip yields a usable active set.

  • drop_reason (str) – Why the trip was dropped, if valid is False.

Parameters:
indices: ndarray
slack: ndarray
margins: ndarray | None
kept_node_rows: ndarray
num_components: int
valid: bool
drop_reason: str = ''
micropurc.core.activeset.active_threshold(x, config)[source]

Compute the flow magnitude below which a link counts as carrying nothing.

The rule is max(threshold_min, threshold_factor * x.max()), so the threshold follows the scale of the flow vector; an empty x gives threshold_min. Route sampling uses this to clean numerically negligible flows before walking the support. Estimation reads activity from dual slack.

Parameters:

  • x (ndarray) – Link flows, shape (L,).

  • config (ActiveSetConfig) – Supplies threshold_factor and threshold_min.

Returns:

The threshold, in flow units.

Return type:

float

micropurc.core.activeset.determine_active_set(x, dest, network, config, projection_mode='robust_pinv', lam_eq=None, c=None, scale_m=None, z_rownorm=None)[source]

Determine the active set for a single trip.

A link is active when \(s_j = c_j + (A^{\top}\lambda)_j\) falls below -config.slack_tol. The active subnetwork is then checked for acyclicity and for reachability of dest, and grounded one node row per connected component; a trip that fails any of these is returned with valid=False and a drop_reason.

Parameters:

  • x (ndarray) – Optimal link flows (unused; kept for call-site symmetry).

  • dest (int) – Destination node id.

  • network (Network) – The network.

  • config (ActiveSetConfig) – Active-set detection settings (slack_tol).

  • projection_mode (str) – Range-space grounding mode (robust_pinv).

  • lam_eq (ndarray | None) – Equality-constraint duals (node potentials) from the solve.

  • c (ndarray | None) – Link cost vector for this trip (Z_variant @ beta).

  • scale_m (ndarray | None) – Per-link perturbation scale. Accepted by the native signature and unread by the slack test.

  • z_rownorm (ndarray | None) – Per-link design-row norms \(\lVert Z_j \rVert\) for this trip’s variant. When given, the returned ActiveSet carries the margins \(|s_j| / \lVert Z_j \rVert\); otherwise margins is None.

Returns:

The trip’s active set.

Raises:

ValueError – If lam_eq, c, or scale_m is missing, or if z_rownorm has a length other than the number of links.

Return type:

ActiveSet

Diagnostics

Convergence and active-set diagnostics for the estimator.

The estimator solves a fixed point whose map is piecewise constant in the active set of binding non-negativity constraints. This module records the per-iteration quantities of a run (residuals, merit values, step sizes, active-set churn, and dual-slack margins), classifies why a run fails to converge, and probes how sensitive the detected active set is to a perturbation of the parameter estimate.

The main entry points are:

class micropurc.core.diagnostics.NonConvergenceReason(value)[source]

Bases: Enum

Classification of why an estimation run failed to converge.

Variables:

  • MAX_ITERATIONS – The iteration budget ran out while the residual was still moving, or the history is too short to classify.

  • STAGNATION – The residual plateaued.

  • OSCILLATION – The residual changes alternate in sign.

  • LINE_SEARCH_EXHAUSTION – The line search rejected most late steps.

  • H_CONDITIONING – The Hessian is ill-conditioned in the late iterations.

  • ACTIVE_SET_CYCLING – Active sets keep flipping while the residual holds still.

  • TRIP_DROP_INSTABILITY – The number of dropped trips fluctuates across iterations.

MAX_ITERATIONS = 'max_iterations'
STAGNATION = 'stagnation'
OSCILLATION = 'oscillation'
LINE_SEARCH_EXHAUSTION = 'ls_exhaustion'
H_CONDITIONING = 'h_conditioning'
ACTIVE_SET_CYCLING = 'aset_cycling'
TRIP_DROP_INSTABILITY = 'trip_drop'
class micropurc.core.diagnostics.IterationRecord(iteration, phi_residual, Q_merit, step_size, accepted, ls_tries, cache_hits=0, cache_misses=0, active_set_changes=0, link_flips=0, cond_H=0.0, min_eig_H=0.0, dropped_trips=0, wall_clock_s=0.0, beta=None, Q_ref=0.0)[source]

Bases: object

Per-iteration diagnostic record.

The residual and merit describe the iterate the estimator entered the iteration with; beta is the iterate it left with.

Variables:

  • iteration (int) – Zero-based iteration counter.

  • phi_residual (float) – \(\lVert \Phi_N(\beta) - \beta \rVert\).

  • Q_merit (float) – Merit value \(Q = \tfrac12 \lVert \Phi_N(\beta) - \beta \rVert^2\).

  • step_size (float) – Damping factor applied to the fixed-point step.

  • accepted (bool) – Whether the line search found a step meeting the non-monotone acceptance test.

  • ls_tries (int) – Number of trial steps evaluated.

  • cache_hits (int) – Filled by callers that track their own caches; the native accumulation path leaves this at 0.

  • cache_misses (int) – Counterpart of cache_hits.

  • active_set_changes (int) – Filled by callers that track per-trip active-set changes; the native accumulation path leaves this at 0.

  • link_flips (int) – Number of trips whose active set differs from the previous iteration, backfilled by DiagnosticsCollector.finalize_churn().

  • cond_H (float) – Condition number of the accumulated Hessian.

  • min_eig_H (float) – Smallest eigenvalue of the accumulated Hessian.

  • dropped_trips (int) – Trips discarded during accumulation at this iterate.

  • wall_clock_s (float) – Seconds spent in this iteration.

  • beta (numpy.ndarray | None) – The iterate after this iteration’s step.

  • Q_ref (float) – Non-monotone reference level the trial step had to beat.

Parameters:
iteration: int
phi_residual: float
Q_merit: float
step_size: float
accepted: bool
ls_tries: int
cache_hits: int = 0
cache_misses: int = 0
active_set_changes: int = 0
cond_H: float = 0.0
min_eig_H: float = 0.0
dropped_trips: int = 0
wall_clock_s: float = 0.0
beta: ndarray | None = None
Q_ref: float = 0.0
class micropurc.core.diagnostics.MarginStats(min_active, median_active, q05_active, q95_active, max_inactive, n_active_total, n_inactive_total)[source]

Bases: object

Summary statistics of dual-slack margins for one iteration.

A link’s distance to flipping is \(|s_j|\). The active statistics are taken over \(|s_j|\) for links with \(s_j < 0\), pooled across the trips of one iteration; small values mark links on the verge of leaving the active set.

Variables:

  • min_active (float) – Smallest \(|s_j|\) over active links, inf when no link is active.

  • median_active (float) – Median \(|s_j|\) over active links.

  • q05_active (float) – 5th percentile of \(|s_j|\) over active links.

  • q95_active (float) – 95th percentile of \(|s_j|\) over active links.

  • max_inactive (float) – Largest slack over inactive links.

  • n_active_total (int) – Number of active links pooled over trips.

  • n_inactive_total (int) – Number of inactive links pooled over trips.

Parameters:
min_active: float
median_active: float
q05_active: float
q95_active: float
max_inactive: float
n_active_total: int
n_inactive_total: int
class micropurc.core.diagnostics.DropReasonCounts(no_active_links=0, trivial_flow=0, dag_connectivity_failed=0, not_dag=0, no_kept_rows=0, other=0)[source]

Bases: object

Per-iteration counts of trip-drop reasons.

Variables:

  • no_active_links (int) – No link passed the dual-slack test.

  • trivial_flow (int) – At most one link is active.

  • dag_connectivity_failed (int) – An active node fails to reach the destination.

  • not_dag (int) – The active subnetwork contains a directed cycle.

  • no_kept_rows (int) – Component grounding left no node row to project on.

  • other (int) – Reasons outside the buckets above.

Parameters:

  • no_active_links (int)

  • trivial_flow (int)

  • dag_connectivity_failed (int)

  • not_dag (int)

  • no_kept_rows (int)

  • other (int)

trivial_flow: int = 0
dag_connectivity_failed: int = 0
not_dag: int = 0
no_kept_rows: int = 0
other: int = 0
property total: int

Total dropped-trip count summed over all reasons.

record(reason)[source]

Increment the counter for a drop reason, bucketing unknowns as other.

Parameters:

reason (str)

Return type:

None

class micropurc.core.diagnostics.DiagnosticsCollector(records=<factory>, active_set_sizes=<factory>, margin_stats=<factory>, drop_reasons=<factory>, peak_rss_kb=0, total_wall_s=0.0, native_time_solve_s=0.0, native_time_projection_s=0.0, native_time_schur_s=0.0, native_total_qp_iters=0, n_threads=0, _prev_trip_active=<factory>, _current_trip_active=<factory>, aset_hash_history=<factory>, per_trip_cumulative_churn=<factory>)[source]

Bases: object

Collects per-iteration diagnostics for a single estimation run.

One collector belongs to one call of fit(), which replaces it at the start of every fit.

Variables:
Parameters:
records: list[IterationRecord]
active_set_sizes: list[list[int]]
margin_stats: list[MarginStats]
drop_reasons: list[DropReasonCounts]
peak_rss_kb: int = 0
total_wall_s: float = 0.0
native_time_solve_s: float = 0.0
native_time_projection_s: float = 0.0
native_time_schur_s: float = 0.0
native_total_qp_iters: int = 0
n_threads: int = 0
aset_hash_history: list[ndarray]
per_trip_cumulative_churn: dict[int, int]
record_iteration(rec)[source]

Append one iteration’s record.

Parameters:

rec (IterationRecord)

Return type:

None

record_active_set_sizes(sizes)[source]

Append the per-trip active-set sizes observed in one iteration.

Parameters:

sizes (list[int])

Return type:

None

record_trip_active_set(trip_id, indices)[source]

Store one trip’s active link indices for the current evaluation.

Parameters:
Return type:

None

record_active_set_hashes(hashes)[source]

Append per-trip active-set hashes from one \(\Phi\) evaluation.

The estimator calls this once before the iteration loop and once per iteration. The arrays are only stored here; finalize_churn() compares them in bulk once the loop is over, which keeps churn accounting off the hot path.

Parameters:

hashes (ndarray) – Per-trip active-set hash, shape (N,). A hash of 0 marks a trip that was dropped.

Return type:

None

flush_active_set_churn(*, detailed=False)[source]

Compare current against previous active sets; return total link flips.

This serves callers that register index sets through record_trip_active_set(). Once any hash array has been recorded, churn is accounted for in bulk by finalize_churn() and this returns zero.

Parameters:

detailed (bool) – Also return the per-trip flip counts.

Returns:

The total number of link flips, or that total paired with a mapping from trip id to flip count when detailed is set.

Return type:

int | tuple[int, dict[int, int]]

finalize_churn()[source]

Compute per-trip and per-iteration churn from the saved hash arrays.

Called once after the iteration loop. Populates per_trip_cumulative_churn and backfills link_flips on each IterationRecord. A trip counts as churning between two evaluations when both of its hashes are nonzero and differ, so a trip dropped at either iterate is passed over.

Return type:

None

record_margins(all_slacks)[source]

Compute and store one iteration’s margin statistics.

Slacks are pooled across trips and split at zero: negative entries are the active links and enter as \(|s_j|\), the rest are the inactive ones.

Parameters:

all_slacks (list[ndarray]) – One dual-slack vector per trip.

Return type:

None

record_drop_reasons(counts)[source]

Append one iteration’s drop-reason counts.

Parameters:

counts (DropReasonCounts)

Return type:

None

snapshot_rss()[source]

Raise the recorded peak RSS to the process’s current peak.

Measurement is best effort: a platform that answers neither getrusage nor the Win32 query leaves the previous value standing.

Return type:

None

peak_rss_mb()[source]

Peak resident set size in MiB, normalized across platforms.

ru_maxrss is bytes on macOS and KiB on Linux, and peak_rss_kb stores ru_maxrss // 1024, so its unit is KiB on macOS and MiB on Linux and Windows. This method converts whichever it holds to MiB.

Returns:

The peak resident set size in MiB.

Return type:

float

accumulate_native_stats(stats)[source]

Sum one native accumulate pass’s phase timings into the fit totals.

The estimator calls the native batch routine once per Krasnoselskii-Mann \(\Phi\) evaluation and once for the covariance pass, and each pass returns its own time_solve_s / time_projection_s / time_schur_s / total_qp_iters. Adding them here makes summary() report totals over the whole fit.

The phase timings are aggregate CPU-seconds summed over the worker threads, so under parallelism they exceed the wall-clock total_wall_s by roughly the achieved thread utilization. They carry the relative cost split across solve, projection, and Schur; total_wall_s carries the runtime.

Parameters:

stats (dict) – One pass’s native statistics dict.

Return type:

None

static compute_cond_H(H)[source]

Compute the condition number and smallest eigenvalue of H.

H is treated as symmetric, and the condition number is \(\lambda_{\max} / \max(|\lambda_{\min}|, 10^{-30})\), the floor keeping a singular H finite.

Parameters:

H (ndarray) – Accumulated Hessian, shape (K, K).

Returns:

The condition number and the smallest eigenvalue, or (inf, 0.0) when the eigenvalue solver fails to converge.

Return type:

tuple[float, float]

classify_non_convergence(stagnation_window=20, stagnation_tol=0.0001, oscillation_window=10, oscillation_frac=0.6, ls_exhaustion_window=10, ls_exhaustion_frac=0.7, cond_threshold=1000000000000.0, churn_window=10, churn_threshold=50, drop_cv_threshold=0.3)[source]

Classify why convergence failed from the iteration history.

The tests run in a fixed order and the first match wins: Hessian conditioning, line-search exhaustion, trip-drop instability, active-set cycling, oscillation, then stagnation. A history that matches none of them, or that holds fewer than three records, is reported as an exhausted iteration budget.

Parameters:

  • stagnation_window (int) – Length of the trailing window shared by the conditioning, line-search, drop, and stagnation tests.

  • stagnation_tol (float) – Relative residual change below which the residual counts as flat, used by the stagnation and cycling tests.

  • oscillation_window (int) – Number of trailing residual differences inspected for sign changes.

  • oscillation_frac (float) – Fraction of those differences that must alternate.

  • ls_exhaustion_window (int) – Shortest trailing window for which the rejection test applies.

  • ls_exhaustion_frac (float) – Fraction of rejected steps that triggers.

  • cond_threshold (float) – Condition number above which the Hessian is blamed.

  • churn_window (int) – Trailing iterations used by the cycling test.

  • churn_threshold (int) – Mean per-iteration link flips that count as churn.

  • drop_cv_threshold (float) – Coefficient of variation of the dropped-trip counts above which drops count as unstable.

Returns:

The reason, paired with an evidence dict holding the metrics that triggered it.

Return type:

tuple[NonConvergenceReason, dict]

summary()[source]

Produce a summary dict for diagnostics output.

The converged entry is a residual-only verdict: the final phi_residual at or below 1e-3. When it comes out False the summary also carries a non_convergence_reason and its evidence. fit() overwrites converged with its own verdict, which uses the configured fixed-point tolerance and the Hessian conditioning bound.

Returns:

The summary, empty when no iteration was recorded.

Return type:

dict

format_convergence_table()[source]

Render the per-iteration convergence diagnostics as a text table.

Returns:

One header line, then one line per iteration carrying the residual, the merit and its non-monotone reference, the step size, the line-search tries and outcome, the cache and churn counters, the dropped trips, the Hessian conditioning, and wall-clock seconds.

Return type:

str

format_activeset_table()[source]

Render the per-iteration active-set diagnostics as a text table.

Rows cover the iterations for which active-set sizes were recorded, so the table is empty when only hashes were collected.

Returns:

One header line, then one line per iteration carrying the trip count, the median, smallest and largest active-set size, the smallest and median active margin, and the churn counter.

Return type:

str

class micropurc.core.diagnostics.ActiveSetReport(converged, num_iterations, final_residual, beta_hat, non_convergence_reason, link_flips_per_iter, residual_per_iter, margin_stats_per_iter, per_trip_cumulative_flips)[source]

Bases: object

Per-run summary of active-set behavior.

Records the per-iteration churn, dual-slack margin distributions, and the non-convergence classification for one estimation run.

Variables:

  • converged (bool) – Whether the run reached the fixed-point tolerance.

  • num_iterations (int) – Total iterations executed.

  • final_residual (float) – \(\lVert \Phi(\beta) - \beta \rVert\) at the last iteration.

  • beta_hat (numpy.ndarray) – The final parameter estimate.

  • non_convergence_reason (str | None) – The NonConvergenceReason value for a run that failed to converge, None otherwise.

  • link_flips_per_iter (list[int]) – Per-iteration total link flips across all trips.

  • residual_per_iter (list[float]) – Per-iteration \(\lVert \Phi(\beta) - \beta \rVert\).

  • margin_stats_per_iter (list[micropurc.core.diagnostics.MarginStats]) – Per-iteration dual-slack margin distributions.

  • per_trip_cumulative_flips (dict[int, int]) – Cumulative link flips per trip, which identifies the trips driving most of the churn.

Parameters:
converged: bool
num_iterations: int
final_residual: float
beta_hat: ndarray
non_convergence_reason: str | None
residual_per_iter: list[float]
margin_stats_per_iter: list[MarginStats]
per_trip_cumulative_flips: dict[int, int]
top_churning_trips(n=20)[source]

Return the n trips with the most cumulative link flips.

Parameters:

n (int)

Return type:

list[tuple[int, int]]

to_dict()[source]

Serialize to a JSON-compatible dict.

Return type:

dict

class micropurc.core.diagnostics.ConvergenceAnalysis(reports)[source]

Bases: object

Aggregate active-set reports across replications and sample sizes.

Cross-tabulates the non-convergence reasons against the sample size and computes summary statistics over replications.

Parameters:

reports (dict[int, list[ActiveSetReport]])

__init__(reports)[source]

Store the per-sample-size replication reports.

Parameters:

reports (dict[int, list[ActiveSetReport]]) – Mapping from sample size N to the list of per-replication reports at that size.

Return type:

None

convergence_rate_by_n()[source]

Fraction of converged runs for each sample size.

Return type:

dict[int, float]

reason_crosstab()[source]

Cross-tabulate the non-convergence reason against sample size.

Returns:

For each N, counts of each non-convergence reason; converged runs are counted under "converged".

Return type:

dict[int, dict[str, int]]

churn_after_iteration(after_iter=3)[source]

Fraction of runs with at least one link flip after a given iteration.

Parameters:

after_iter (int) – Only count flips from this iteration onward.

Returns:

For each N, the fraction of runs with a link flip after after_iter.

Return type:

dict[int, float]

bias_comparison(beta_true)[source]

Compare estimator bias between converged and non-converged runs.

Parameters:

beta_true (ndarray) – The true parameter vector.

Returns:

For each N, a dict with the converged and non-converged bias and RMSE arrays.

Return type:

dict[int, dict[str, ndarray]]

format_reason_table()[source]

Render the non-convergence reason by sample-size cross-tab as a table.

Return type:

str

class micropurc.core.diagnostics.StabilityProbeResult(n_trips, n_directions, radii, trips_changed, links_flipped)[source]

Bases: object

Result of the active-set stability probe.

Variables:

  • n_trips (int) – Number of trips probed.

  • n_directions (int) – Random directions drawn per radius.

  • radii (list[float]) – Perturbation radii, in the units of the parameter norm.

  • trips_changed (dict[float, list[int]]) – Radius to a per-trip list holding the largest active-set symmetric difference seen over the directions; a direction that makes the trip invalid contributes 1.

  • links_flipped (dict[float, list[int]]) – Radius to the symmetric-difference sizes pooled over all trips and directions that changed at all.

Parameters:
n_trips: int
n_directions: int
radii: list[float]
trips_changed: dict[float, list[int]]
format_table()[source]

Format the probe result as a text table over perturbation radii.

Return type:

str

micropurc.core.diagnostics.active_set_stability_probe(network, solver, beta_hat, y, b, active_set_config, projection_mode='robust_pinv', radii=None, n_directions=5, seed=42)[source]

Perturb \(\hat\beta\) and measure how stable the active set is.

For every radius and every random unit direction the estimate is perturbed, the forward problem is re-solved for each trip, and the resulting active set is compared with the one at \(\hat\beta\). Trips whose reference active set is invalid are skipped, and a perturbed solve that turns invalid is scored as a change of size one.

Parameters:

  • network – Network supplying the incidence and the design Z.

  • solver – Forward solver for the inner QP.

  • beta_hat (ndarray) – Parameter estimate to perturb, shape (K,).

  • y (ndarray) – Route indicators, shape (N, L); only the trip count is read.

  • b (ndarray) – Node imbalances, shape (N, V); each destination is its argmax.

  • active_set_config – Active-set detection settings.

  • projection_mode (str) – Range-space grounding mode.

  • radii (list[float] | None) – Perturbation radii. Defaults to [1e-4, 1e-3, 1e-2, 1e-1] scaled by the norm of beta_hat.

  • n_directions (int) – Random unit directions drawn per radius.

  • seed (int) – Seed for drawing the directions.

Returns:

The per-radius active-set changes and flip counts.

Return type:

StabilityProbeResult

Lazy design

Lazy, term-coded design matrix for large problems.

A LazyDesign holds per-link attribute columns, per-variant covariates, and a compact table of terms that the native engine evaluates on the fly, so the (N, L, K) design is never held in memory. Each term contributes one parameter’s column, evaluated per link and per trip variant:

  • TERM_COLUMN(a) — link attribute column a;

  • TERM_PRODUCT(a, b) — product of columns a and b (interaction);

  • TERM_COVARIATE_INDICATOR(a, cov, op, thr) — column a gated by the indicator 1[covariate[cov] op thr].

Link attributes may vary across an attribute-variant axis (for instance, time-of-day); variant_attr maps each trip variant to its attribute-variant row, and variant_covariates holds the covariate values used by indicator terms. Build one directly, or via micropurc.spec.compile_design() with materialize=False.

class micropurc.core.lazy_design.LazyDesign(link_attrs, variant_attr, variant_covariates, trip_variant, term_code, term_a, term_b, term_cov, term_op, term_thr)[source]

Bases: object

A term-coded design evaluated lazily by the native engine.

Variables:

  • link_attrs (numpy.ndarray) – Attribute columns, shape (n_attr_variants, L, A).

  • variant_attr (numpy.ndarray) – Attribute-variant row for each variant, shape (n_var,).

  • variant_covariates (numpy.ndarray) – Covariate values per variant, shape (n_var, C).

  • trip_variant (numpy.ndarray) – Variant index for each trip, shape (N,).

  • term_code (numpy.ndarray) – Opcode per term, shape (K,).

  • term_a (numpy.ndarray) – Primary attribute column per term, shape (K,).

  • term_b (numpy.ndarray) – Second column for product terms, shape (K,).

  • term_cov (numpy.ndarray) – Covariate index for indicator terms, shape (K,).

  • term_op (numpy.ndarray) – Comparison code for indicator terms, shape (K,).

  • term_thr (numpy.ndarray) – Threshold for indicator terms, shape (K,).

Parameters:
variant_attr: ndarray
variant_covariates: ndarray
trip_variant: ndarray
term_code: ndarray
term_a: ndarray
term_b: ndarray
term_cov: ndarray
term_op: ndarray
term_thr: ndarray
property num_attributes: int

Number of parameters (one per term).

Number of links.

property num_variants: int

Number of trip variants.