Copyright (c) 2026 - Pierre M. Fiorini, Ph.D. Busy-period grammar computation and paper-fragment generator. This module is the executable companion to the paper "On Computing the Mean Stationary Response Time in M/M/1 Fork-Join $(n,k)$ Queues" It provides two deliberately separate command-line workflows. Copyright (c) 2026 - Pierre M. Fiorini, Ph.D. ===============================================================================TWO EXECUTION OPTIONS=============================================================================== Option 1: regular computation, no LaTeX generation-------------------------------------------------- python busy_period_grammar_fresh.py compute \ --n 3 --k 2 --lambda-rate 0.4 --mu 1 1 1 \ --tolerance 1e-10 --max-expansions 500000 This mode runs the analytic best-first frontier computation for one model andprints the usual numerical diagnostics. It does not write tables, figures,PGFPlots code, or paper fragments. Add --des to run an independent-replicationDES validation for the same model. This is the normal one-case computationalworkflow. Option 2: paper generation with fresh values-------------------------------------------- python busy_period_grammar_fresh.py paper \ --out generated_fresh \ --des-reps-small 30 --des-warmup-small 5000 --des-measured-small 50000 \ --des-reps-large 20 --des-warmup-large 3000 --des-measured-large 20000 This mode recomputes all analytic table entries from the best-first probabilitymass frontier solver, runs fresh DES replications, and then writes every table,TikZ diagram, and PGFPlots graph fragment used in the paper. Unlike earlierversions of this generator, paper mode does not populate validation tables fromfrozen manuscript constants. The displayed analytic values, DES means, DEShalf-widths, absolute differences, pending masses, and expansion counts are allwritten from the current execution. For quick debugging, paper mode accepts --skip-des, --case-filter, and small DESreplication settings. For manuscript-grade regeneration, use the DES settingsreported in the paper or larger values. ===============================================================================MODEL IMPLEMENTED BY THE ANALYTIC FRONTIER SOLVER=============================================================================== The analytic model is the infinite-buffer Markovian cancel-on-k fork-join queue.Jobs arrive as a Poisson process of rate lambda. Each job creates one task ateach of n FCFS M/M/1 branches. Branch i has exponential service rate mu_i. Thejob departs when any k of its n copies complete service, and all unfinishedcopies of that same job are cancelled immediately. The state used by the analytic solver is a queue-length state, not a job-idstate. This is valid because all arrivals are synchronous and all branches areFCFS. A horizontal queue depth represents one active job. If a servicecompletion occurs at branch i and the branch length before completion is r, thenumber of unfinished copies of that depth-r job is M_r(q) = #{j : q_j >= r}. (eq:Mr) After the mechanical completion, cancellation fires exactly when M_r(q) - 1 = r}. (eq:cancelvec) When cancellation does not fire, only the completing branch is decremented.These formulas are implemented by full_transitions() and grouped_transitions(). The full heterogeneous state is q=(q_1,...,q_n). The total rate out of anon-empty state is Sigma(q) = lambda + sum_{i:q_i>0} mu_i. (eq:sigma) The active-job count is the number of occupied horizontal queue depths: N_sys(q) = max_i q_i. (eq:Nsysfull) For homogeneous or partially homogeneous systems, branches are partitioned intorate classes. For rate class g, c_{g,a} is the number of branches in group gwith queue length a. The grouped state is x = (c_1,...,c_G), sum_a c_{g,a} = n_g. (eq:groupstate) The grouped exit rate and active-job count are Sigma(x) = lambda + sum_g mu_g sum_{a>=1} c_{g,a}, (eq:groupsigma) N_sys(x) = max{a : sum_g c_{g,a} > 0}. (eq:groupsigma) The grouped rank count is M_a(x) = sum_h sum_{b>=a} c_{h,b}. (eq:groupMr) Grouping is exact symmetry lumping; it is not a mean-field approximation and itis not a finite-buffer truncation. ===============================================================================RENEWAL-REWARD EQUATIONS AND FRONTIER EXPANSION=============================================================================== The empty state is regenerative. Starting immediately after the initiatingarrival, let tau_0 be the first return time to the empty state. The steady-statemean job response time is computed by the renewal-reward identity E[T] = E[ integral_0^{tau_0} N_sys(X(t)) dt ] / E[A(tau_0)]. (eq:rr-intro) The first-step residual reward equations are R(x) = N_sys(x)/Sigma(x) + sum_{y!=0} P(x,y) R(y), (eq:Rfirst) J(x) = lambda/Sigma(x) + sum_{y!=0} P(x,y) J(y). (eq:Jfirst) The initiating arrival is counted separately. Therefore the displayed responseratio is E[T] = R(1) / (1 + J(1)). (eq:grammarrr) AnalyticalBestFirstFrontierSolver evaluates these equations forward over thebusy-period grammar. It stores a sparse pending frontier measure m(x). Poppinga state with mass m contributes m * N_sys(x)/Sigma(x) to accumulated job-time reward, contributes m * lambda/Sigma(x) to expected future admitted jobs, and pushes each nonterminal successor with itsembedded-chain probability. Empty successors are absorbed. The traversal order is best-first: the largest pending mass is expanded next.This is an evaluation order, not a different stochastic model. If the frontieris fully drained with no pruning and no resampling, the result equals the exactanalytic regenerative ratio in exact arithmetic and is numerically exact up tofloating-point roundoff. ===============================================================================ERROR-ACCOUNTING EQUATIONS=============================================================================== After finitely many expansions, the current non-empty frontier F has residuals R_F = sum_{x in F} m(x) R(x), J_F = sum_{x in F} m(x) J(x). (eq:residuals) With no pruning or resampling, T = (R_hat + R_F) / (J_hat + J_F). (eq:frontier_exact) Subtracting the prefix estimate T_hat=R_hat/J_hat gives T - T_hat = (R_F J_hat - R_hat J_F)/(J_hat(J_hat+J_F)). (eq:error_decomp) Pending mass alone is not an error bound on an unbounded state space. If U_Rand U_J are nonnegative reward envelopes satisfying the paper's inequalities,then the computable residual interval is R_hat/(J_hat + Jbar_F) <= T <= (R_hat + Rbar_F)/J_hat. (eq:interval) When deterministic frontier resampling is used, the moved mass is D_move = 1/2 sum_x |F_tilde(x) - F(x)|, (eq:moved) and bounded test functions satisfy the total-variation perturbation bound |sum_x f(x)F_tilde(x)-sum_x f(x)F(x)| <= osc_F(f) D_move. (eq:tvbound) The mass-conservation diagnostic is epsilon_mass = |m_absorbed + m_pending + m_pruned - 1|. (eq:masserr) ===============================================================================DES VALIDATION AND CONFIDENCE INTERVALS=============================================================================== The script includes an event-level DES simulator with the same FCFS cancel-on-ksemantics as the analytic grammar. The simulator uses job identifiers, branchqueues, per-branch in-service tokens, and a global event heap. Tokens invalidatestale service-completion events whenever an in-service residual copy is cancelledby a job completion on another branch. For R independent replications, after warm-up, replication r has mean Tbar_r = (1/M) sum_{j=1}^M T_{rj}. (eq:repmean) The reported DES confidence interval is Tbar_DES +/- t_{0.975,R-1} s_R/sqrt(R), s_R^2 = (1/(R-1)) sum_r (Tbar_r - Tbar_DES)^2. (eq:ci) The DES interval is a simulation sampling interval only. The analytic grammarvalue is deterministic and is accompanied by frontier diagnostics rather than asampling interval.