api.fj

fj_respt_nosplit(K, lambda, mu)

Mean response time of the distributed no-splitting parallel system

In the distributed no-splitting policy a job is a sequence of K tasks that is routed in one piece to a single server chosen uniformly among the K servers, so each server sees a Poisson stream of rate lambda/K whose service time is the sum of K exponential stages of rate mu. That is an M/E_K/1 queue, and the Pollaczek-Khinchine mean response time reduces to

R(rho) = 1/(mu - lambda), rho = lambda/mu, R_{D/NS} = [ K - (K-1)*rho/2 ] * R(rho),

which is the reference against which the splitting policies of the same section are judged: splitting wins because it leaves fewer servers idle. At K = 1 the expression collapses to the M/M/1 response time.

Parameters:

K – Number of servers, equal to the number of tasks per job
lambda – Total job arrival rate
mu – Task service rate at each server, mu > lambda for stability

Returns:

R – Mean job response time under distributed no splitting rho: Utilization of each server, lambda/mu

Examples

R = fj_respt_nosplit(K, lambda, mu)
[R, rho] = fj_respt_nosplit(K, lambda, mu)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Section 6.1 on page 17:33.

Original: R. Nelson, D. Towsley, A. N. Tantawi, “Performance Analysis of Parallel Processing Systems”, IEEE Trans. Software Eng. 14(4), 1988.

fj_respt_closed(K, x, M, A)

Varki bound on the residence time of a closed fork-join subnetwork

Bounds the mean residence time of a K-way fork-join subnetwork of homogeneous exponential servers with mean service time x. Writing A for the mean number of jobs an arriving job finds at the subnetwork,

R_{P_K}(M) <= x * [ H_K + A ].

In a closed network the arrival theorem supplies A = Q(M-1); when the network consists of the parallel subsystem alone, every one of the other M-1 jobs is necessarily inside it, so A = M-1 and

R_{P_K}(M) <= x * [ H_K + M - 1 ],

which holds with equality for K = 2 and is therefore exact there.

Parameters:

K – Number of parallel servers (positive integer)
x – Mean service time of each server
M – Number of jobs circulating in the closed network
A – Mean queue length seen on arrival (optional, default M-1)

Returns:

R – Upper bound on the mean residence time of the subnetwork exact: True when the bound is known to be tight, i.e. K = 2 and A = M-1

Examples

R = fj_respt_closed(K, x, M)
[R, exact] = fj_respt_closed(K, x, M, A)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Eq. (67) and Theorem 4.1 on page 17:47.

Original: E. Varki, “Response Time Analysis of Parallel Computer and Storage Systems”, IEEE Trans. Parallel Distrib. Syst. 12(11), 2001.

fj_respt_bulk(K, lambda, mu, c, nmax)

Centralized splitting analysed as an M[K]/M/c bulk arrival system

Under centralized splitting a request forks into K tasks that are held in a single central queue and served by c identical servers, so the same server may serve several tasks of the same request. That is an M[K]/M/c queue with fixed batch size K. Its stationary distribution has no product form, and is obtained here by solving the balance equations of the level chain truncated at nmax, which is exact up to the tail mass that truncation discards.

Two response times are reported. The task response time follows from Little’s law on the mean number in system. The request response time is the completion of the LAST of the K tasks of a tagged request: by PASTA the batch finds n tasks in system, its last task is the (n+K)-th in line, and under first come first served with c exponential servers it starts service after max(0, n+K-c) departures, each an exponential of rate c*mu, whence

E[R_request] = sum_n p_n * [ max(0, n+K-c)/(c*mu) + 1/mu ].

This is the centralized counterpart of the distributed splitting fork-join system, and lower bounds it because no task is bound to a particular server.

Parameters:

K – Batch size, that is the number of tasks per request
lambda – Arrival rate of requests (batches)
mu – Task service rate at each server
c – Number of servers
nmax – Truncation level of the task-count chain (optional)

Returns:

Rreq – Mean request response time, that is the last of the K tasks Rtask: Mean task response time Q: Mean number of tasks in system p: Stationary distribution of the number of tasks in system

Examples

[Rreq, Rtask] = fj_respt_bulk(K, lambda, mu, c)
[Rreq, Rtask, Q, p] = fj_respt_bulk(K, lambda, mu, c, nmax)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Section 6.1 on page 17:33.

Original: R. Nelson, D. Towsley, A. N. Tantawi, “Performance Analysis of Parallel Processing Systems”, IEEE Trans. Software Eng. 14(4), 1988.

fj_qgb(D, P, M, Z)

Geometric bound on the queue length of a fork-join subnetwork

Non-iterative geometric bound on the mean queue length of each fork-join subnetwork of a closed queueing network. Subnetwork n consists of P(n) parallel queues traversed as a P(n)-way fork-join request, and carries a per-visit service demand D(n).

y_n(M) = D_n * M / (Z + sum_j D_j * H_{P_j} + Dmax * M) Q_n(M) = H_{P_n} * [ y_n/(1-y_n) - y_n^(M+1)/(1-y_n) ]

with Dmax = max_j D_j and H_k the k-th harmonic number. The harmonic weights are what distinguishes this from the ordinary geometric bound of pfqn_qzgblow: a P-way fork-join subnetwork inflates its own demand by H_P in the denominator and its queue length by H_P in the numerator. Setting P(n) = 1 for every n recovers pfqn_qzgblow exactly.

Parameters:

D – Vector of per-visit service demands, one per subnetwork
P – Vector of fork degrees, one per subnetwork (P(n) >= 1)
M – Number of circulating jobs (positive integer)
Z – Think time (optional, default 0)

Returns:

Q – Vector of bounded mean queue lengths, one per subnetwork y: Vector of geometric ratios y_n(M)

Examples

[Q, y] = fj_qgb(D, P, M)
[Q, y] = fj_qgb(D, P, M, Z)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Eq. (70).

Original: G. Casale, R. R. Muntz, G. Serazzi, “Geometric Bounds: A Noniterative Analysis Technique for Closed Queueing Networks”, IEEE Trans. Computers 57(6), 2008.

fj_ordstat_exp(ri, k)

M = FJ_ORDSTAT_EXP(RI, K)

Mean of the K-th smallest of N independent EXPONENTIAL branch completion times whose means are RI, i.e. the instant a k-of-n (quorum) join fires. K = N is the ordinary AND-join, the maximum, and K = 1 the minimum.

With lambda_i = 1/RI(i) and m = N-K stragglers allowed,

E[X_(K)] = sum_{j=m+1..N} (-1)^(j-m-1) C(j-1,m) e_j, e_j = sum_{|S|=j} 1 / sum_{i in S} lambda_i

the inclusion-exclusion identity for the order statistics of independent exponentials. At m = 0 it collapses to sum_j (-1)^(j-1) e_j, the classical expression for the maximum, TERM BY TERM: a full join therefore evaluates exactly as it did before this function existed.

The sum has 2^N terms and its signs alternate, so it is evaluated exactly only while the branch count is small. Beyond MAXEXACT branches a genuine quorum (K < N) is evaluated by FJ_QUORUM_MOMENTS instead, whose Poisson-binomial recurrence adds no cancellation; a full join keeps the exact path at every N so that no existing result moves.

Reference: A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys 47(2), Article 17, 2014, Sec. 3 (Eq. 18-19).

fj_lst_max_het(lambda, s)

Laplace-Stieltjes transform of the maximum of heterogeneous exponentials

Evaluates L*_K(s) = E[exp(-s Y)] for Y = max(X_1,…,X_K) with independent X_i ~ Exp(lambda_i), using the Harrison and Zertal recurrence over the sub-collections obtained by deleting one rate at a time:

( s + sum_{j=1..m} lambda_j ) L*_m(lambda, s) = sum_{j=1..m} lambda_j * L*_{m-1}(lambda j, s), 1 <= m <= K,

anchored at L*_0 = 1 because the maximum of an empty collection is zero. The recurrence is evaluated bottom-up over the 2^K sub-collections, each identified by a bit mask, so every value is computed once.

Parameters:

lambda – Vector of K positive exponential rates
s – Transform argument, scalar or vector, s >= 0

Returns:

L – Transform value, same shape as s

Examples

L = fj_lst_max_het(lambda, s)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Eq. (29).

Original: P. G. Harrison, S. Zertal, “Queueing Models of RAID Systems with Maxima of Waiting Times”, Performance Evaluation 64(7-8), 2007.

fj_ism_green(lambda, mu, s, c)

Green’s independent server model of simultaneous server requests

In the independent server model a customer requires j servers simultaneously to begin service, with probability c(j), and then releases them asynchronously as each of its j tasks completes at rate mu. Servers can sit idle while customers wait, which is what separates the model from M/G/s, and customer service ends when the last of its tasks does, so the customer service time is the maximum of j exponentials, H_j/mu.

The cycle decomposition of Green splits time into a queueing period, during which a queue exists, and a nonqueue period. Everything follows from two quantities:

E[B] = sum_j c(j) * sum_{i=0..j-1} 1/((s-i)*mu),

the interservice time, which is the j-th order statistic of s exponentials because all s servers are busy whenever a customer enters service during a queueing period, and

E[D] = sum_i sum_{k=1..i} [ sum_{m=0..k-1} 1/((i-m)*mu) ] * q(i)*c(s-i+k)/p_d,

the initial delay of the customer that starts a queueing period, which finds i servers busy and must wait for k of them to free. The busy-server distribution q during a nonqueue period and its mean length come from the embedded chain of arrivals and completions absorbed when a queue forms:

V = (I - T)^-1, E[Qbar] = sum_j v(s,j)/(lambda + j*mu), q(i) = v(s,i)/((lambda + i*mu) * E[Qbar]), E[Q] = E[D]/(1 - lambda*E[B]), p_q = E[Q]/(E[Q] + E[Qbar]).

The waiting time transform of Eq. (61) factors into the equilibrium transform of D and the Pollaczek-Khinchine transform of an M/G/1 queue with service B, so its mean is

E[W] = (1 - pi0) * [ E[D^2]/(2*E[D]) + lambda*E[B^2]/(2*(1-lambda*E[B])) ], pi0 = (1 - lambda*E[B]) / (1 - lambda*(E[B] - E[D])).

The second moments of B and D are exact, each stage of the order statistic being an independent exponential.

Note that Eq. (65) of the survey prints the inner sum of E[D] as starting at 1/(s*mu) even though only i servers are busy; the sum is started at 1/(i*mu) here, which is what the accompanying text prescribes and what makes E[D] reduce to E[B] when i = s.

Parameters:

lambda – Customer arrival rate
mu – Per-task service rate
s – Number of servers
c – Vector of s probabilities, c(j) = P(a customer needs j servers)

Returns:

W – Mean waiting time before service starts R: Mean response time, waiting plus the maximum of the j tasks out: Struct with EB, EB2, ED, ED2, EQ, EQbar, pq, pd, pi0, rho, q, ES

Examples

[W, R, out] = fj_ism_green(lambda, mu, s, c)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Section 6.4, Eqs. (59)-(65).

Original: L. Green, “A Queueing System in Which Customers Require a Random Number of Servers”, Operations Research 28(6), 1980.

fj_is_homogeneous(sn)

Checks for a single fork-join pair with homogeneous parallel branches

This is NOT a test for the presence of Fork/Join nodes (use sn_has_fork_join for that): most fork-join models fail this predicate. It tests membership in the homogeneous class on which the FJ_codes tail approximation of Qiu, Perez and Harrison (IFIP Performance 2015) is defined, i.e. the requirements: - Single Fork-Join pair - K parallel queues between Fork and Join - Homogeneous service distributions across parallel queues - Supported distributions (Exp, HyperExp(2), Erlang(2), MAP(2)) - Open classes only

Parameters:: sn – Network structure
Returns:: isHomogeneous – True if network is valid FJ topology for FJ_codes fjInfo: Struct with fields: forkIdx, joinIdx, queueIdx, K, errorMsg

Examples

[isHomogeneous, fjInfo] = fj_is_homogeneous(sn)

References

fj_dispersion(shape, rate, d, options)

Mean subtask dispersion of a split-merge system with Erlang branches

Subtask dispersion is the interval between the first and the last subtask completion of a request. For heterogeneous branches with distribution functions F_i shifted by deterministic delays d_i, the heterogeneous order statistics of Eq. (24) give

E[X_(N)] = integral_0^inf [ 1 - prod_i F_i(x - d_i) ] dx, E[X_(1)] = integral_0^inf prod_i [ 1 - F_i(x - d_i) ] dx, E[D_d] = E[X_(N)] - E[X_(1)] = integral_0^inf [ 1 - prod_i F_i(x-d_i) - prod_i (1-F_i(x-d_i)) ] dx.

The integrand of the last line is what is evaluated; it is non-negative and vanishes at both ends, whereas the difference of the two products printed in the survey is not the dispersion and can go negative.

Branch i is an Erlang with shape(i) stages of rate rate(i), which is the split-merge equivalent used in the delay-scheduling construction: a subtask with q other subtasks ahead of it in its parallel queue behaves as an Erlang(q+1, mu). Quadrature is composite Simpson on a horizon widened until the completion probability is within TOL of one.

Parameters:

shape – Vector of N Erlang stage counts (positive integers)
rate – Vector of N Erlang stage rates (positive)
d – Vector of N non-negative deterministic delays (optional, default zeros)
options – Optional struct with fields tol and npoints

Returns:

Edisp – Mean subtask dispersion Emax: Mean completion time of the last subtask Emin: Mean completion time of the first subtask

Examples

Edisp = fj_dispersion(shape, rate)
[Edisp, Emax, Emin] = fj_dispersion(shape, rate, d)
[Edisp, Emax, Emin] = fj_dispersion(shape, rate, d, options)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Eqs. (24)-(25) and Section 6.2.

Original: I. Tsimashenka, W. J. Knottenbelt, “Reduction of Subtask Dispersion in Fork-Join Systems”, EPEW, 2013.

fj_delay_opt(shape, rate, options)

Deterministic subtask delays that minimise mean dispersion

Chooses the vector of deterministic delays d = (d_1,…,d_N) that minimises the mean subtask dispersion E[D_d] of fj_dispersion. Holding back a fast branch costs little at the last completion and buys a great deal at the first, so the minimiser is generally interior and strictly positive on every branch but the slowest.

The objective is minimised by cyclic coordinate descent with a golden section line search on each coordinate, which is deterministic and needs no derivative; the delay of the branch with the largest mean is pinned at zero, because adding a constant to every delay shifts both order statistics equally and leaves the dispersion unchanged. The search therefore returns the representative with min(d) = 0.

Parameters:

shape – Vector of N Erlang stage counts (positive integers)
rate – Vector of N Erlang stage rates (positive)
options – Optional struct with fields tol, npoints, maxsweeps and dtol

Returns:

d – Vector of N optimal non-negative delays with min(d) = 0 Edisp: Mean dispersion attained at d Emax: Mean completion time of the last subtask at d, the price paid

Examples

d = fj_delay_opt(shape, rate)
[d, Edisp, Emax] = fj_delay_opt(shape, rate, options)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Section 6.2.

Original: I. Tsimashenka, W. J. Knottenbelt, “Reduction of Subtask Dispersion in Fork-Join Systems”, EPEW, 2013.

fj_dag_makespan(pred, rate)

Makespan of a task system with precedence constraints

A task system is a set of n tasks whose precedence relations form a directed acyclic graph. Because the precedence relation is acyclic, so is the chain whose state is the SET of completed tasks, and the makespan can be swept level by level instead of solved as a linear system.

In a state whose completed set is S, the tasks eligible to run are those all of whose predecessors lie in S; a task i among the k of them completes at rate rate(i,k), so the state is held for M(S) = 1/T(S) with T(S) the sum of those rates, and moves to S union {i} with probability b(S,i) = rate(i,k)/T(S). Making the rate depend on how many tasks run concurrently is what couples the task system to the queueing network underneath it: two tasks sharing a processor each run slower than either would alone.

The weighted delay to reach a state and the probability of reaching it obey

p(R) = sum_{S -> R} p(S)*b(S,R), D(R) = M(R)*p(R) + sum_{S -> R} b(S,R)*D(S),

started at p(empty) = 1, and the makespan is D at the fully completed state. The initiation, completion and execution times of each task follow from the same sweep, D being accumulated over the transitions that start and that finish the task.

Parameters:

pred – n by n logical matrix, pred(i,j) true when task i precedes task j
rate – Completion rates: an n-vector of constant rates, or an n by n matrix whose entry (i,k) is the rate of task i while k tasks run concurrently

Returns:

C – Mean makespan, the time until every task has completed I: Vector of mean initiation times, one per task Cend: Vector of mean completion times, one per task E: Vector of mean execution times, Cend - I

Examples

C = fj_dag_makespan(pred, rate)
[C, I, Cend, E] = fj_dag_makespan(pred, rate)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Eq. (66) and Section 7.4 on page 17:49.

Original: A. Thomasian, P. F. Bay, “Analytic Queueing Network Models for Parallel Processing of Task Systems”, IEEE Trans. Computers C-35(12), 1986.

fj_cox_fit(m1, c2)

Two-stage Coxian fit of a mean and a squared coefficient of variation

Marie’s balanced-stage fit of a two-stage Coxian law to a target mean 1/mu and squared coefficient of variation c2. The two stages are required to contribute equally to the mean, 1/mu1 = q/mu2, which closes the system of two moment equations in three unknowns and yields

mu1 = 2*mu, q = 1/(2*c2), mu2 = 2*mu*q = mu/c2.

The fit needs q <= 1, hence c2 >= 0.5; below that the balanced-stage condition is infeasible and an Erlang stage count is the natural choice instead. The admissible number of exponential stages of an Erlang or same-rate Coxian representation of the same target is bracketed by

ceil(1/c2) <= k <= floor(1/c2) + 1,

which is also returned.

Parameters:

m1 – Target mean (positive)
c2 – Target squared coefficient of variation, c2 >= 0.5

Returns:

mu1 – Rate of the first stage mu2: Rate of the second stage q: Probability that the second stage is visited kmin: Smallest admissible Erlang stage count for the same c2 kmax: Largest admissible Erlang stage count for the same c2

Examples

[mu1, mu2, q] = fj_cox_fit(m1, c2)
[mu1, mu2, q, kmin, kmax] = fj_cox_fit(m1, c2)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Section 4.5, Eqs. (39)-(40) and the discussion of Algorithm 3.2.5.

Original: R. Marie, “Methodes iteratives de resolution de modeles mathematiques de systemes informatiques”, RAIRO Informatique 12(2), 1978; A. O. Allen, “Probability, Statistics and Queueing Theory”, Academic Press, 2nd ed., 1990, Algorithm 3.2.5.

fj_char_max_discrete(K, dist_type, param)

Characteristic maximum of a discrete random variable

Gravey’s characteristic maximum for lattice distributions. Let m_K be the smallest integer with P(X > m_K) <= 1/K; then

M_K = m_K + K * sum_{k >= m_K} P(X > k),

which upper bounds the expected maximum of K i.i.d. copies of X and costs O(1) instead of the alternating binomial sum. Two lattice laws admit a closed form for the tail sum:

geometric, P(X = k) = (1-p) p^k: m_K = ceil( -ln(K) / ln(p) ), M_K = m_K + K * p^(m_K+1) / (1-p), and the exact maximum is E[Y_K] = sum_{k=1..K} binom(K,k) (-1)^(k+1) p^k / (1 - p^k);

Poisson, P(X = k) = exp(-theta) theta^k / k!: M_K = m_K * (1 - K*P(X > m_K)) + K*theta*P(X > m_K - 1), which is the same tail sum rewritten through E[(X-m)^+] = theta*P(X>m-1) - m*P(X>m).

The exact expected maximum is returned alongside: in closed form for the geometric law, and by summing 1 - F(k)^K over the lattice for the Poisson.

Parameters:

K – Number of i.i.d. copies (positive integer)
dist_type – ‘geometric’ or ‘poisson’
param – Success-complement p in (0,1) for the geometric, mean theta > 0 for the Poisson

Returns:

MK – Characteristic maximum, an upper bound on the expected maximum mK: Smallest integer with P(X > mK) <= 1/K exact: Exact expected maximum of the K copies

Examples

[MK, mK] = fj_char_max_discrete(K, dist_type, param)
[MK, mK, exact] = fj_char_max_discrete(K, dist_type, param)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Eq. (48) and the lattice cases on page 17:25.

Original: A. Gravey, “A Simple Construction of an Upper Bound for the Mean of the Maximum of N Identically Distributed Random Variables”, J. Applied Probability 22(4), 1985.

fj_char_max_blom(K, Finv, alpha, beta)

Blom-corrected plotting position for the characteristic maximum

The characteristic maximum m_K is the quantile at which the survival function of a branch drops to 1/K. The naive plotting position m_K = F^-1(K/(K+1)) is biased; Blom’s correction replaces it with

m_K = F^-1( (K - alpha) / (K - alpha - beta + 1) ),

which for alpha = beta = 0 falls back on the naive position. The survey quotes alpha = 0.4886 and beta = 0.3140, which are the defaults here.

For the standard normal branch the position is bracketed without any inversion, for K >= 5, by

sqrt(2 ln K - ln ln K - 3) < m_K < sqrt(2 ln K - ln ln K),

and the leading term alone gives the Kruskal-Weiss estimate m_K ~ mu + sigma sqrt(2 ln K).

Parameters:

K – Number of i.i.d. copies (positive integer)
Finv – Quantile function handle (optional; the standard normal by default)
alpha – Blom numerator offset (optional, default 0.4886)
beta – Blom denominator offset (optional, default 0.3140)

Returns:

mK – Blom-corrected characteristic maximum lo: Kruskal-Weiss lower bracket, valid for the standard normal and K >= 5 hi: Kruskal-Weiss upper bracket, valid for the standard normal and K >= 5

Examples

mK = fj_char_max_blom(K)
[mK, lo, hi] = fj_char_max_blom(K)
[mK, lo, hi] = fj_char_max_blom(K, Finv, alpha, beta)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Section 4.9 on page 17:26.

Original: G. Blom, “Statistical Estimates and Transformed Beta Variables”, Wiley, 1958; W. Kruskal, G. Weiss, “Allocation of Observations in Comparison of Treatments”, 1985.

sn_fj_validate(sn)

SN_FJ_VALIDATE(SN)

Validate that a fork-join model is within the supported feature set of the native CTMC/SSA fork-join implementation (v1): closed classes only, non-nested fork-join pairs, standard join strategy (wait for all siblings), one task per output link. Unsupported models raise a line_error with a specific message; further structural checks (class switching on a branch, nesting, sibling traps) are performed during the branch discovery in ModelAdapter.fjtag.

fj_xmax_erlang(K, k, mu)

Expected maximum of K i.i.d. Erlang-k random variables

Computes the expected value of the maximum of K independent and identically distributed Erlang random variables with k stages and rate mu per stage.

For k=2 (Erlang-2): X_K^max = (1/mu) * sum_{n=1}^{K} C(K,n) * (-1)^{n-1} * sum_{m=1}^{n} C(n,m) * m! / (2*n^{m+1})

For general k-stage Erlang, the formula involves k-fold convolution of exponential distributions.

The Erlang distribution has mean k/mu and variance k/mu^2, giving CV = 1/sqrt(k) < 1.

Parameters:

K – Number of parallel servers (positive integer)
k – Number of Erlang stages (positive integer, default=2)
mu – Rate parameter per stage (mean service time = k/mu)

Returns:

Xmax – Expected maximum service time X_K^max

Examples

Xmax = fj_xmax_erlang(K, k, mu)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Page 17:13.

fj_xmax_emma(K, param, dist_type)

Expected maximum using EMMA method

Computes the expected value of the maximum of K i.i.d. random variables using the EMMA (Expected Maximum from Marginal Approximation) method.

The method is based on the property that: [F_X(E[Y_K])]^K ≈ φ = 0.570376

Therefore: E[Y_K] ≈ F^{-1}(φ^{1/K})

For the exponential distribution with rate μ: E[Y_K] = -(1/μ) * ln(1 - φ^{1/K})

The constant φ = exp(-exp(-γ)) ≈ 0.570376, where γ ≈ 0.5772 is the Euler-Mascheroni constant.

Parameters:

K – Number of random variables (positive integer)
param – Rate mu for exponential, or inverse CDF function handle
dist_type – Optional: ‘exp’ (default) or ‘general’

Returns:

Xmax – Approximate expected maximum E[Y_K]

Examples

Xmax = fj_xmax_emma(K, mu)                    % Exponential with rate mu
Xmax = fj_xmax_emma(K, F_inv)                 % General CDF inverse
Xmax = fj_xmax_emma(K, mu, 'exp')             % Explicit exponential
Xmax = fj_xmax_emma(K, F_inv, 'general')      % Explicit general

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Eq. (44) on page 17:22.

Original: Y. Sun, K.L. Peterson, “Computing Extreme Order Statistics from Large Data Sets”, 2012.

ge_fit(ET, VT)

Two-moment fit of a generalized exponential law

Matches the generalized exponential CDF

F(x) = (1 - exp(-x/beta))^alpha

on a mean and a variance,

E[T] = beta*(psi(alpha+1) - psi(1)), V[T] = beta^2*(psi’(1) - psi’(alpha+1)).

The squared coefficient of variation depends on the SHAPE ALONE and decreases monotonically in it, so the shape is recovered by a scalar root-find on a logarithmic scale and the scale then follows in closed form. SCV = 1 is the exponential case alpha = 1, kept exact.

It is the per-branch law of the ForkTail approximation, shared by FJ_TAIL_FORKTAIL (the AND-join) and FJ_TAIL_ORDSTAT (the k-of-n quorum), which is why it lives in a file of its own rather than inside either.

Parameters:

ET – Mean of the law, positive
VT – Variance of the law, positive

Returns:

alpha – Fitted shape parameter beta: Fitted scale parameter

Examples

[alpha, beta] = ge_fit(ET, VT)

References

M. Nguyen, S. Alesawi, N. Li, H. Che, H. Jiang, “ForkTail: A Black-Box Fork-Join Tail Latency Prediction Model for User-Facing Datacenter Workloads”, ACM HPDC 2018, pp. 206-217.

fj_rmax_evd(K, R, sigma_R, calibrated)

Maximum response time using Extreme Value Distribution approximation

Approximates the expected value of the maximum response time for K parallel servers using the Extreme Value Distribution (EVD).

R_K^max(rho) = R(rho) + (sqrt(6)*ln(K)/pi) * sigma_R(rho)

where R(rho) is the mean response time and sigma_R(rho) is the standard deviation of response time.

This approximation is based on fitting the response time distribution to a Gumbel (Type I extreme value) distribution with: - Location parameter: a = R - gamma*b - Scale parameter: b = sqrt(6)*sigma_R/pi

where gamma ≈ 0.5772 is the Euler-Mascheroni constant.

Note: Thomasian et al. [2007] found that dividing the correction term by 1.27 improves accuracy.

Parameters:

K – Number of parallel servers (positive integer)
R – Mean response time
sigma_R – Standard deviation of response time
calibrated – Optional: use calibrated formula (default: false)

Returns:

Rmax – Approximate maximum response time R_K^max

Examples

Rmax = fj_rmax_evd(K, R, sigma_R)
Rmax = fj_rmax_evd(K, R, sigma_R, calibrated)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Eq. (43) on page 17:21.

fj_xmax_moments_het(lambda, n)

Moments of the maximum of heterogeneous exponentials by recurrence

Exact moments of Y = max(X_1,…,X_K) with independent X_i ~ Exp(lambda_i), obtained by differentiating the Harrison and Zertal transform recurrence n times and evaluating at the origin:

M_m(lambda, n) = [ n * M_m(lambda, n-1) + sum_{j=1..m} lambda_j * M_{m-1}(lambda j, n) ] / sum_{j=1..m} lambda_j,

with M_m(lambda, 0) = 1 and M_0(., n) = 0 for n >= 1. Every moment up to order n is produced by the same sweep, so all of them are returned.

Eq. (30) of the survey prints the second sum WITHOUT the lambda_j weight. That form is not the derivative of Eq. (29) and does not reproduce the textbook two-variable answer 1/l1 + 1/l2 - 1/(l1+l2); the weight is restored here, and fj_xmax_het gives the independent inclusion-exclusion check.

Parameters:

lambda – Vector of K positive exponential rates
n – Highest moment order (optional, default 1)

Returns:

M – Row vector of moments of orders 1..n of the maximum

Examples

M = fj_xmax_moments_het(lambda)
M = fj_xmax_moments_het(lambda, n)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Eq. (30).

Original: P. G. Harrison, S. Zertal, “Queueing Models of RAID Systems with Maxima of Waiting Times”, Performance Evaluation 64(7-8), 2007.

fj_rmax_erlang(K, k, lambda, mu)

Maximum response time R_K^max for Erlang service times

Computes the expected value of the maximum response time for K M/E_k/1 queues (Poisson arrivals, Erlang-k service times).

For K=2 (two servers), uses the closed-form formula: R_2^max = R_1 + R_2 - sum_{m=0}^{k1-1} sum_{n=0}^{k2-1} C(m+n,m) * mu1^m * mu2^n / (mu1+mu2)^{m+n+1}

For general K, uses numerical integration: R_K^max = integral_0^inf [1 - prod_{i=1}^{K} F_Erlang(t)] dt

where F_Erlang(t) = 1 - exp(-mu*t) * sum_{j=0}^{k-1} (mu*t)^j / j!

Parameters:

K – Number of parallel servers (positive integer)
k – Number of Erlang stages (positive integer)
lambda – Arrival rate
mu – Service rate per Erlang stage (mean service = k/mu)

Returns:

Rmax – Expected maximum response time R_K^max

Examples

Rmax = fj_rmax_erlang(K, k, lambda, mu)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Eq. (33) and (34) on page 17:18.

fj_rmax(K, lambda, mu)

Compute maximum response time R_K^max(rho) for K M/M/1 queues

Computes the expected value of the maximum of K response times from independent M/M/1 queues. This serves as an upper bound to the K-way Fork-Join response time R_K^{F/J}(rho).

R_K^max(rho) = H_K * R(rho) = H_K / (mu - lambda)

where H_K is the K-th Harmonic number and R(rho) is the M/M/1 mean response time.

This formula holds because for M/M/1 queues, the response time distribution is exponential with rate (mu - lambda), and the expected maximum of K i.i.d. exponential random variables with rate theta is H_K / theta.

Parameters:

K – Number of parallel servers (positive integer)
lambda – Arrival rate
mu – Service rate (mu > lambda for stability)

Returns:

Rmax – Expected maximum response time R_K^max(rho) = H_K * R(rho)

Examples

Rmax = fj_rmax(K, lambda, mu)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Page 17:8, Eq. (1), and surrounding text.

fj_xmax_hz_het(m1, m2, cdf, options)

Harrison-Zertal approximation of the maximum of general variables

Expected maximum of K independent but not identically distributed non-negative variables, each supplied through its first two moments and its distribution function. Writing S for a sub-collection and alpha_i = 1/m1_i, the recurrence averages, over which branch is singled out, the expected maximum of the remaining branches plus the residual life still owed by the branch singled out:

I(S) = (1/|S|) * sum_{i in S} [ I(S i) + (m2_i / (2*m1_i)) * L*_{Si}(alpha_i) ],

anchored at I({i}) = m1_i. The transform of the maximum over a sub-collection is recovered from the product of the distribution functions,

L*_T(s) = s * integral_0^inf exp(-s t) * prod_{j in T} F_j(t) dt,

evaluated by composite Simpson quadrature on a horizon widened until the product of the distribution functions is within TOL of one.

For identically distributed branches the recurrence collapses onto the closed form of fj_xmax_hz, and for identical exponential branches it is exact at H_K/lambda.

Parameters:

m1 – Vector of K branch means (positive)
m2 – Vector of K branch second moments, m2(i) >= m1(i)^2
cdf – Cell array of K function handles, cdf{i}(t) = P(X_i <= t)
options – Optional struct with fields tol and npoints

Returns:

Xmax – Approximate expected maximum

Examples

Xmax = fj_xmax_hz_het(m1, m2, cdf)
Xmax = fj_xmax_hz_het(m1, m2, cdf, options)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Eq. (46).

Original: P. G. Harrison, S. Zertal, “Queueing Models of RAID Systems with Maxima of Waiting Times”, Performance Evaluation 64(7-8), 2007.

fj_respt_vm(K, lambda, mu)

Varma-Makowski approximation for K-way F/J response time

Computes an approximate mean response time for a K-way Fork-Join queueing system with Poisson arrivals and exponential service times.

R_K^{F/J}(rho) ≈ [H_K + (A_K - H_K) * rho] * (mu - lambda)^{-1}

where A_K = sum_{i=1}^{K} C(K,i) * (-1)^{i-1} * sum_{m=1}^{i} C(i,m) * (m-1)! / i^{m+1}

and C(n,k) denotes the binomial coefficient “n choose k”.

Parameters:

K – Number of parallel servers (positive integer)
lambda – Arrival rate
mu – Service rate (mu > lambda for stability)

Returns:

R – Approximate K-way F/J response time R_K^{F/J}(rho)

Examples

R = fj_respt_vm(K, lambda, mu)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Table I, Eq. (4) on page 17:10.

Original: S. Varma and A.M. Makowski, “Interpolation Approximations for Symmetric Fork-Join Queues”, Performance Evaluation, 20, 1994.

fj_xmax_hz(m1, m2, K)

Harrison-Zertal approximation of the maximum of i.i.d. variables

Closed form obtained by collapsing the Harrison and Zertal recurrence onto K identically distributed branches described by their first two moments:

X_K^max ~ m1 + ( m2 / (2*m1) ) * ( H_K - 1 ).

The correction is the equilibrium (residual life) mean of the branch distribution scaled by H_K - 1, so it reads as “one branch, plus the expected residual work still owed by the branches that finish later”. Writing m2/(2*m1) = m1*(1+SCV)/2 gives the equivalent m1 * [ 1 + (1+SCV)/2 * (H_K - 1) ], which is exact for the exponential distribution (SCV = 1) and reduces to m1 at K = 1 for every branch law.

Parameters:

m1 – Mean of the branch distribution (positive)
m2 – Second moment of the branch distribution, m2 >= m1^2
K – Number of branches (positive integer)

Returns:

Xmax – Approximate expected maximum resid: Equilibrium mean m2/(2*m1) used as the per-branch inflation

Examples

Xmax = fj_xmax_hz(m1, m2, K)
[Xmax, resid] = fj_xmax_hz(m1, m2, K)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Eq. (46) and the identically distributed specialization on page 17:23.

Original: P. G. Harrison, S. Zertal, “Queueing Models of RAID Systems with Maxima of Waiting Times”, Performance Evaluation 64(7-8), 2007.

fj_respt_varki(K, lambda, mu)

Varki et al. approximation for K-way F/J response time

Computes an approximate mean response time for a K-way Fork-Join queueing system with Poisson arrivals and exponential service times. This is the mean of pessimistic (upper) and optimistic (lower) bounds.

R_K^{F/J}(rho) ≈ (1/mu) * [H_K + (rho/(2*(1-rho))) * (S1 + (1-2*rho)*S2)]

where: S1 = sum_{i=1}^{K} 1/(i - rho) S2 = sum_{i=1}^{K} 1/(i*(i - rho))

Parameters:

K – Number of parallel servers (positive integer)
lambda – Arrival rate
mu – Service rate (mu > lambda for stability)

Returns:

R – Approximate K-way F/J response time R_K^{F/J}(rho)

Examples

R = fj_respt_varki(K, lambda, mu)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Table I, Eq. (5) on page 17:10.

Original: E. Varki, A. Merchant, J. Xu, and X. Qiu, “Issues and Challenges in the Performance Analysis of Real Disk Arrays”, IEEE Trans. Parallel Distrib. Syst., 25(6), 2014 (published as Varki et al. [2013] in draft).

fj_xmax_het(lambda, n)

Exact moments of the maximum of heterogeneous exponentials

Exact n-th moment of Y = max(X_1,…,X_K) for independent but not identically distributed exponential variables, X_i ~ Exp(lambda_i), by inclusion-exclusion on the survival function:

E[Y^n] = sum over the nonempty subsets S of {1..K} of (-1)^(|S|+1) * n! / ( sum_{i in S} lambda_i )^n .

For n = 1 and K = 2 this collapses to the textbook 1/l1 + 1/l2 - 1/(l1+l2), and for equal rates to H_K/lambda.

The cost is 2^K - 1 terms, so the enumeration is refused beyond K = 24; use fj_xmax_moments_het for the O(2^K) but numerically safer recursion, or fj_xmax_exp when the rates are equal.

Parameters:

lambda – Vector of K positive exponential rates
n – Moment order (optional, default 1)

Returns:

Mn – n-th moment of the maximum

Examples

Xmax = fj_xmax_het(lambda)
Mn = fj_xmax_het(lambda, n)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Eq. (27) and the inclusion-exclusion expansion of Section 7.3 on page 17:47.

fj_respt_nt(K, lambda, mu)

Nelson-Tantawi approximation for K-way F/J response time

Computes an approximate mean response time for a K-way Fork-Join queueing system with Poisson arrivals and exponential service times. Valid for 2 <= K <= 32.

R_K^{F/J}(rho) ≈ [H_K/H_2 + (1 - H_K/H_2) * 4*rho/11] * (1.5 - rho/8) / (mu - lambda)

The approximation is based on the exact 2-way solution and a scaling approximation that upper and lower bounds increase at the same rate. The coefficient alpha(rho) ≈ 4*rho/11 was obtained from simulation.

Parameters:

K – Number of parallel servers (2 <= K <= 32)
lambda – Arrival rate
mu – Service rate (mu > lambda for stability)

Returns:

R – Approximate K-way F/J response time R_K^{F/J}(rho)

Examples

R = fj_respt_nt(K, lambda, mu)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Table I, Eq. (3) on page 17:10.

Original: R. Nelson and A.N. Tantawi, “Approximate Analysis of Fork/Join Synchronization in Parallel Queues”, IEEE Trans. Computers, 37(6), 1988.

fj_respt_2way(lambda, mu)

Exact two-way Fork-Join response time R_2^{F/J}(rho)

Computes the exact mean response time for a 2-way (K=2) Fork-Join queueing system with Poisson arrivals and exponential service times. This is derived from the analysis in Flatto and Hahn [1984].

R_2^{F/J}(rho) = (H_2 - rho/8) * R(rho) = (12 - rho)/8 * R(rho)

where H_2 = 1.5 and R(rho) = (mu - lambda)^{-1} is the M/M/1 mean response time.

Note that R_2^{F/J}(rho) < R_2^max(rho), with the difference being rho/8 * R(rho).

Parameters:

lambda – Arrival rate
mu – Service rate (mu > lambda for stability)

Returns:

R – Exact 2-way F/J response time R_2^{F/J}(rho)

Examples

R = fj_respt_2way(lambda, mu)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Eq. (6) on page 17:10.

Original derivation: L. Flatto and S. Hahn, “Two Parallel Queues Created by Arrivals with Two Demands I”, SIAM J. Appl. Math., 44(5), 1984.

fj_quantile(K, q, F_inv)

Quantile approximation for maximum of K random variables

Approximates the q-th quantile of the maximum of K i.i.d. random variables from the standard Gumbel (Type I extreme value) distribution.

x(K,q) ≈ ln(K) - ln(ln(1/q))

where P(M_K <= x(K,q)) = q is the q-th quantile.

Note: This approximation is inaccurate for small values of K.

For general distributions, the quantile of the maximum can be computed using the inverse CDF: x(K,q) = F^{-1}(q^{1/K})

Parameters:

K – Number of random variables (positive integer)
q – Quantile probability (0 < q < 1)
F_inv – Optional: inverse CDF function handle for general distributions

Returns:

x_Kq – q-th quantile of the maximum of K random variables

Examples

x_Kq = fj_quantile(K, q)                  % Standard Gumbel approx
x_Kq = fj_quantile(K, q, F_inv)           % General distribution

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Page 17:22.

fj_order_stat(y, k, K, F_X)

CDF and expected value of k-th order statistic

Computes the CDF and expected value of the k-th order statistic (k-th smallest) of K i.i.d. random variables.

For maximum (k=K): F_{Y_K}(y) = [F_X(y)]^K

For k-th order statistic (k-th smallest): F_{Y_k}(y) = sum_{j=k}^{K} C(K,j) * F_X(y)^j * (1-F_X(y))^{K-j}

Expected value of maximum: E[Y_K] = integral_0^inf [1 - F_X(y)^K] dy

Parameters:

y – Value(s) at which to evaluate CDF
k – Order of the statistic (1 = minimum, K = maximum)
K – Total number of random variables
F_X – CDF function handle: F_X(y) returns CDF value

Returns:

F_Yk – CDF of k-th order statistic at y E_Yk: Expected value of k-th order statistic (if requested)

Examples

F_Yk = fj_order_stat(y, k, K, F_X)
[F_Yk, E_Yk] = fj_order_stat(y, k, K, F_X)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Eq. (18) and (19) on page 17:15.

fj_harmonic(K)

Compute Harmonic sum H_K = sum(1/k) for k=1 to K

Computes the K-th Harmonic number, which is the sum of reciprocals from 1 to K. This is a fundamental quantity in Fork-Join analysis.

For large K, H_K ≈ ln(K) + γ, where γ ≈ 0.57721 is the Euler-Mascheroni constant.

Parameters:: K – Number of parallel servers (positive integer)
Returns:: H – Harmonic sum H_K = 1 + 1/2 + 1/3 + … + 1/K

Examples

H = fj_harmonic(K)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014.

fj_gk_bound(K, type)

Compute G(K) factors for expected maximum approximation

Computes the G(K) factor used in the approximation: X_K^max ≈ mu_X + sigma_X * G(K)

For standardized distributions (mean=0, variance=1), G(K) = X_K^max.

Available G(K) formulas: - Exponential: G(K) = H_K - 1 - Uniform: G(K) = sqrt(3) * (K-1) / (K+1) - EVD: G(K) = sqrt(6) * ln(K) / pi - Upper bound: G(K) <= (K-1) / sqrt(2K-1) (tight bound from David [1970])

Parameters:

K – Number of random variables (positive integer)
type – Optional: ‘exp’, ‘uniform’, ‘evd’, ‘bound’, or ‘all’ (default)

Returns:

GK – G(K) value(s) - scalar or struct depending on type

Examples

GK = fj_gk_bound(K)                    % Returns struct with all
GK = fj_gk_bound(K, 'exp')             % Exponential G(K)
GK = fj_gk_bound(K, 'uniform')         % Uniform G(K)
GK = fj_gk_bound(K, 'evd')             % EVD G(K)
GK = fj_gk_bound(K, 'bound')           % Upper bound

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Eq. (22) and (23) on page 17:16.

Original: H.A. David, “Order Statistics”, Wiley, 1970.

fj_mg1_respt_moments(lambda, ES, ES2, ES3)

Mean and variance of the M/G/1 response time, as ForkTail inputs

Closes the white-box route of fj_tail_forktail for a fork branch that is an M/G/1 FCFS queue, from the first three moments of its service time:

E[T] = E[S]*(1 + rho/(1-rho) * (1+SCV_S)/2) V[T] = E[W]^2 + lambda*E[S^3]/(3*(1-rho)) + E[S^2] - E[S]^2

with rho = lambda*E[S] and E[W] = lambda*E[S^2]/(2*(1-rho)), i.e. the Pollaczek-Khinchine mean waiting time. The third moment enters only the variance, so an exponential or phase-type branch needs no extra input beyond what the service distribution already reports (getMoments(3) on a Markovian law, or 3*ES*VS + ES^3 + skewness*VS^(3/2) in general).

Parameters:

lambda – Arrival rate at the branch
ES – First moment of the service time
ES2 – Second moment of the service time
ES3 – Third moment of the service time

Returns:

ET – Mean response time VT: Variance of the response time

Examples

[ET, VT] = fj_mg1_respt_moments(lambda, ES, ES2, ES3)

References

M. Nguyen, S. Alesawi, N. Li, H. Che, H. Jiang, “ForkTail: A Black-Box Fork-Join Tail Latency Prediction Model for User-Facing Datacenter Workloads”, ACM HPDC 2018, Eqs. (10) and (11).

fj_char_max(K, param, dist_type)

Characteristic maximum M_K for order statistics

Computes the characteristic maximum M_K, which provides a bound for the expected value of the maximum of K i.i.d. random variables.

Let m_K be the greatest lower bound such that P(X > m_K) <= 1/K. For continuous distributions with a density: P(X > m_K) = 1/K.

The characteristic maximum is: M_K = m_K + K * integral_{m_K}^{inf} P(X > x) dx

For specific distributions: - Exponential: m_K = ln(K)/μ, M_K = H_K/μ - Erlang-k: m_K solves exp(-μ*m_K) * sum_{i=0}^{k-1} (μ*m_K)^i/i! = 1/K M_K = (k/μ) * [1 + K*exp(-μ*m_K)*(μ*m_K)^k/k!]

Parameters:

K – Number of random variables (positive integer)
param – Distribution parameters (rate mu, or [k, mu] for Erlang)
dist_type – Optional: ‘exp’ (default), ‘erlang’, or ‘general’

Returns:

MK – Characteristic maximum M_K mK: Threshold m_K where P(X > m_K) = 1/K

Examples

[MK, mK] = fj_char_max(K, mu)                     % Exponential
[MK, mK] = fj_char_max(K, mu, 'exp')              % Exponential (explicit)
[MK, mK] = fj_char_max(K, [k, mu], 'erlang')      % Erlang-k
[MK, mK] = fj_char_max(K, S_X, 'general')         % General (S_X is survival function)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Eq. (47) on page 17:24.

A. Gravey, “A Simple Construction of an Upper Bound for the Mean of the Maximum of n Identically Distributed Random Variables”, 1985.

fj_bounds(K, lambda, mu)

Upper and lower bounds for K-way F/J response time

Computes pessimistic (upper) and optimistic (lower) bounds for the mean response time of a K-way Fork-Join queueing system.

Upper bound (pessimistic): R_K^max(rho) = H_K / (mu * (1 - rho)) = H_K * R(rho)

Lower bound (optimistic) from Varki et al.: R_K^{F/J(opt)}(rho) = (1/mu) * [H_K + S_{K(K-rho)}] where S_{K(K-rho)} = sum_{j=1}^{K} (1/j) * (rho/(j - rho))

The bounds satisfy: R_K^{F/J(opt)}(rho) <= R_K^{F/J}(rho) <= R_K^max(rho)

Parameters:

K – Number of parallel servers (positive integer)
lambda – Arrival rate
mu – Service rate (mu > lambda for stability)

Returns:

Rmax – Upper bound (pessimistic): R_K^max(rho) Rmin: Lower bound (optimistic): R_K^{F/J(opt)}(rho)

Examples

[Rmax, Rmin] = fj_bounds(K, lambda, mu)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Eq. (1) and Eq. (2) on page 17:9.

fj_tsm_capacity(s, f, r, x)

Saturation throughput of the team service model

In the team service model a class-k job seizes r(k) of the s servers at once, holds them for a mean x(k), and releases them all together. The apparent saturation rate is the one that drives the server utilization to one,

Lambda_max = s / sum_k f(k)*r(k)*x(k),

but it is attainable only when the scheduler can pack jobs into execution states that leave no server idle. The attainable capacity is the largest arrival rate for which some mixture p over the feasible execution states balances every class,

maximise Lambda subject to sum_j p_j * n(j,k) / x(k) = Lambda * f(k), sum_j p_j = 1, p >= 0,

over the execution states j, each a multiset of jobs with total server demand at most s. n(j,k) counts the class-k jobs in state j. The linear program is solved by a dense simplex on its standard form; the optimum equals Lambda_max exactly when every state carrying positive probability is full capacity.

For the two-server two-class special case with r = (1,2), the strict first come first served discipline cannot pack at all and reaches only

lambda_FCFS = 2*mu1*mu2 / (f1^2*mu2 + 2*f2^2*mu1 + 2*f1*f2*(mu1+mu2)),

which is returned whenever the arguments describe that case.

Parameters:

s – Number of servers
f – Vector of class frequencies in the arrival stream, summing to one
r – Vector of per-class server requirements (positive integers, r(k) <= s)
x – Vector of per-class mean service times (positive)

Returns:

Lmax – Apparent saturation rate that drives the utilization to one Llp: Attainable capacity from the linear program, Llp <= Lmax Lfcfs: Strict first come first served capacity, NaN outside the two-server two-class case p: Execution-state probabilities attaining Llp, as a struct with fields states and prob

Examples

Lmax = fj_tsm_capacity(s, f, r, x)
[Lmax, Llp, Lfcfs, p] = fj_tsm_capacity(s, f, r, x)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Section 6.3, Eqs. (57)-(58) and the linear programming formulation on page 17:36.

Original: A. Thomasian, “Performance Evaluation of Centralized Databases with Static Locking”, IEEE Trans. Software Eng. SE-11(4), 1985; K. Omahen, L. Schrage, “A Queueing Analysis of a Multiprocessor System with Shared Memory”, Symp. Computer-Communication Networks, 1972.

fj_simplex(A, b, c)

Two-phase dense simplex for small equality-form linear programs

Solves min c’z subject to A z = b and z >= 0 with a dense two-phase simplex using Bland’s rule, which terminates without cycling at the cost of a slower pivot sequence. The programs raised by the fork-join capacity analyses have at most a few hundred columns, so the dense tableau is the right trade and keeps the four language ports free of any linear programming dependency.

Phase one drives a full set of artificial variables out of the basis; if their sum cannot be driven to zero the program is infeasible. Phase two then optimises the true objective on the remaining basis.

Parameters:

A – Equality constraint matrix, m by n
b – Right hand side, m by 1
c – Cost row, 1 by n, minimised

Returns:

z – Optimal point, 1 by n status: 0 optimal, 1 infeasible, 2 unbounded, 3 iteration limit obj: Optimal objective value

Examples

[z, status] = fj_simplex(A, b, c)
[z, status, obj] = fj_simplex(A, b, c)

sn_fj_foldback(QN, UN, RN, TN, CN, XN, fjclassmap, Korig)

[QN,UN,RN,TN,CN,XN] = SN_FJ_FOLDBACK(QN,UN,RN,TN,CN,XN,FJCLASSMAP,KORIG)

Fold the auxiliary-class columns of the average metrics computed on an FJ tag-augmented struct (ModelAdapter.fjtag) back into the original classes: queue lengths, utilizations and throughputs of the sibling classes are exact aggregates of the original class they were forked from; response times are recomputed by Little’s law after folding.

QN,UN,RN,TN are (nstations x Kaug); CN,XN are (1 x Kaug); the outputs retain only the first Korig columns.

sn_build_fj_sync_map(sn)

FJSYNCMAP = SN_BUILD_FJ_SYNC_MAP(SN)

Builds a fork-join synchronization map from LINE’s sn structure. For each (Fork, Join) pair, identifies which source nodes feed into the join and must be synchronized using mmap_max.

Output:

fjSyncMap.nodeSync(joinIdx, srcIdx) = groupId

groupId > 0 means srcIdx belongs to sync group groupId at joinIdx
groupId == 0 means srcIdx is an independent (non-synced) flow

fjSyncMap.forkOfGroup(groupId) = forkIdx fjSyncMap.joinOfGroup(groupId) = joinIdx fjSyncMap.nGroups = total number of sync groups

fj_quorum_moments(branchMeans, branchVars, k)

FJ_QUORUM_MOMENTS Mean and variance of a k-of-n (quorum) join completion time.

[M, V] = FJ_QUORUM_MOMENTS(BRANCHMEANS, BRANCHVARS, K) returns the first two moments of the K-th smallest of N independent branch completion times, each branch being given by its mean and variance.

Each branch is expanded into a discrete step CDF by a two-moment fit. At each time the number of completed branches is Poisson-binomial, so its distribution is built by the recurrence q_j <- q_{j-1}*F_i + q_j*(1-F_i) over branches i, and the K-th order statistic is the upper tail sum_{j>=K} q_j. For K = N this reduces to the product of the branch CDFs, i.e. the ordinary AND-join, and for K = 1 to 1 - prod(1 - F_i), i.e. the minimum.

This is the same quantity as the inclusion-exclusion identity: F_(k)(t) = sum_{i=k..n} (-1)^(i-k) * C(i-1, k-1) * e_i(F_1(t), …, F_n(t))

used by LQNS, but every term here is a probability in [0,1] and none is subtracted, so it avoids the catastrophic cancellation the alternating binomial sum incurs as N grows.

Follows the formulation of Omari, Franks, Woodside and Pan, as implemented in LQNS 6.x (randomvar.cc). Mirrors jline.api.fj.FJ_quorum in the Java runtime.

fj_serialization(Rs, R0, M, alpha)

Blocking probability and pseudoserver delay of serialization phases

A serialization phase is a stretch of a job’s execution that is protected by an exclusive lock, so at most one of the M circulating jobs may occupy it. The queueing network is no longer product form, and the delay in entering phase s is represented by a pseudoserver that is bypassed when the phase is free. A job entering phase s is blocked when at least one of the other M-1 jobs is inside it, and treating those jobs as independently placed in proportion to the residence times gives

P_s(M) = 1 - [ 1 - R_s(M)/R(M) ]^(M-1), R(M) = sum_s R_s(M),

where the sum runs over the nonserialized phase and every serialization phase. The delay charged at the pseudoserver is alpha*R_s(M), with alpha depending on the residence time distribution and on where within the phase the blocked job arrives; alpha = 1/2 is the value for an arrival uniform in a phase of low utilization, which is the regime in which the approximation is stated.

Parameters:

Rs – Vector of mean residence times inside each serialization phase
R0 – Mean residence time in the nonserialized phase
M – Number of circulating jobs (positive integer)
alpha – Fraction of the phase charged to a blocked job (optional, default 0.5)

Returns:

P – Vector of probabilities of being blocked on entering each phase delay: Vector of expected pseudoserver delays, P(s)*alpha*Rs(s) Rtot: Mean cycle time including the serialization delays

Examples

P = fj_serialization(Rs, R0, M)
[P, delay, Rtot] = fj_serialization(Rs, R0, M, alpha)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Section 7.1 on page 17:43.

Original: A. Thomasian, “Queueing Network Models to Estimate Serialization Delays in Computer Systems”, Performance, 1983.

fj_xmax_pareto(K, beta, k)

Expected maximum for Pareto distribution

Computes the expected value of the maximum of K i.i.d. Pareto random variables using order statistics formulas.

Pareto distribution CDF: F(x) = 1 - α(x + γ)^{-β}, x >= 0, β > 2

For the standardized Pareto with mean 1: α = γ^β, γ = β - 1 Mean = 1, Second moment M = 2 + 2/(β - 2)

The j-th moment of the i-th order statistic is: m_{i,j} = Γ(K+1) * Γ(K - j + 1 - i/α) / [Γ(K + 1 - i/α) * Γ(K - j + 1)] * k^i

For the maximum (i = K), the expected value can be computed numerically or approximated using the characteristic maximum method.

Parameters:

K – Number of random variables (positive integer)
beta – Shape parameter (beta > 2 required for finite moments)
k – Optional: scale parameter (default: beta - 1 for mean = 1)

Returns:

Xmax – Expected maximum E[Y_K] MK: Characteristic maximum M_K

Examples

Xmax = fj_xmax_pareto(K, beta)               % Standard Pareto (mean=1)
Xmax = fj_xmax_pareto(K, beta, k)            % Pareto with scale k
[Xmax, MK] = fj_xmax_pareto(K, beta, k)      % Also return characteristic max

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Page 17:24 and 17:26.

N.L. Johnson, S. Kotz, N. Balakrishnan, “Continuous Univariate Distributions”, Vol. 1, Wiley, 1995.

fj_xmax_2(lambda1, lambda2)

Expected maximum of 2 exponential random variables

Computes the expected value of the maximum of 2 independent exponential random variables with rates lambda1 and lambda2.

Y_2^max = 1/lambda1 + 1/lambda2 - 1/(lambda1 + lambda2)

This is the exact closed-form solution for K=2 case.

For identical rates (lambda1 = lambda2 = lambda): Y_2^max = 2/lambda - 1/(2*lambda) = 3/(2*lambda) = 1.5/lambda = H_2/lambda

which matches the general formula X_K^max = H_K/mu for exponentials.

Parameters:

lambda1 – Rate of first exponential
lambda2 – Rate of second exponential (optional, default=lambda1)

Returns:

Xmax – Expected maximum Y_2^max

Examples

Xmax = fj_xmax_2(lambda1, lambda2)
Xmax = fj_xmax_2(lambda)  % Same rate for both

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Eq. (27) on page 17:17.

fj_synch_delay(lambda, mu)

Synchronization delay S_2(rho) for two-way Fork-Join

Computes the mean synchronization delay for a 2-way Fork-Join queueing system. This is the average time tasks spend waiting at synchronization queues for their siblings to complete.

S_2(rho) = (1/2) * (1 - rho/4) * R(rho)

The F/J response time decomposes as: R_2^{F/J}(rho) = R(rho) + S_2(rho)

where R(rho) is the delay at the server (M/M/1 response time).

Parameters:

lambda – Arrival rate
mu – Service rate (mu > lambda for stability)

Returns:

S – Mean synchronization delay S_2(rho)

Examples

S = fj_synch_delay(lambda, mu)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Eq. (8) on page 17:11.

sn_fj_supports(sn)

[BOOL, REASON] = SN_FJ_SUPPORTS(SN)

@brief Can the exact fork-join construction be asked for this model?

The fork-join model class SN_FJ_VALIDATE admits, asked as a predicate rather than raised. SolverCTMC.supportsModelMethod and SolverSSA.supportsModelMethod call it so that a CALLER (model.help, findSolver, SolverAUTO) sees the verdict before paying for a run, and BOTH analyzers reach the SAME rules through MODELADAPTER.FJTAG -> SN_FJ_VALIDATE. The sentence the validator raises names ‘the native CTMC/SSA fork-join implementation’, which is why this predicate lives beside it rather than inside either solver.

IT WRAPS THE VALIDATOR RATHER THAN RESTATING IT, and that is the point: the rules are eight and they move (pairing, join strategy, tasks-per-link, branch probability, open classes through a fork), so a second copy here would be a second thing to keep in step. There is exactly one body of rules and two ways in – one that raises, for the run, and this one, which answers.

WHAT IT REFUSES AND WHY THE ANALYZER IS RIGHT TO. The fork-join PAIRING is a declaration carried by the Join (joinOf here, _fork in native python, the third constructor argument in all four codebases), not a derivation from the routing: a nested model such as examples/basic/forkJoin/fj_basic_nesting has two forks and two joins whose pairing the routing alone does not determine. So a Join built without naming its fork leaves SN.FJ empty, and ‘Fork nodes without a matched Join’ is the honest answer to a model that declares none – not a topology test that failed to see one.

@param sn NetworkStruct of the model @return bool true when the fork-join construction may run @return reason the refusal, or ‘’ when BOOL is true

fj_xmax_normal(K, mu, sigma, method)

Expected maximum for normal distribution

Computes the expected value of the maximum of K i.i.d. normal random variables using the approximation from Johnson et al. [1995].

E[Y_K] ≈ μ + σ * [sqrt(2*ln(K)) - (ln(ln(K)) - ln(4π) + 2γ) / (2*sqrt(2*ln(K)))]

where γ ≈ 0.5772 is the Euler-Mascheroni constant.

A bias correction from Petzold [2000] can be applied: δ(K) = 0.1727 * K^(-0.2750)

The variance of the maximum is approximately: Var[Y_K] ≈ 1.64492 * σ^2 / (2*ln(K))

A simpler approximation from Arnold [1980] is also available: E[Y_K] ≈ μ + σ * sqrt(2*ln(K))

Parameters:

K – Number of random variables (positive integer >= 2)
mu – Mean of the normal distribution
sigma – Standard deviation of the normal distribution
method – Optional: ‘johnson’ (default), ‘arnold’, or ‘corrected’

Returns:

Xmax – Approximate expected maximum E[Y_K] Vmax: Approximate variance of maximum Var[Y_K]

Examples

Xmax = fj_xmax_normal(K, mu, sigma)
Xmax = fj_xmax_normal(K, mu, sigma, 'johnson')     % Johnson et al. (default)
Xmax = fj_xmax_normal(K, mu, sigma, 'arnold')      % Arnold approximation
Xmax = fj_xmax_normal(K, mu, sigma, 'corrected')   % With Petzold correction
[Xmax, Vmax] = fj_xmax_normal(...)                 % Also return variance

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Eq. (45) on page 17:23.

N.L. Johnson, S. Kotz, N. Balakrishnan, “Continuous Univariate Distributions”, Vol. 2, Wiley, 1995.

fj_extract_params(sn, fjInfo)

Extract FJ_codes parameters from LINE network structure

Extracts arrival process, service process, and K value from a validated Fork-Join network structure for use with FJ_codes analysis.

Parameters:

sn – Network structure (after sn_nonmarkov_toph conversion)
fjInfo – FJ topology info from fj_is_homogeneous (forkIdx, joinIdx, queueIdx, K)

Returns:

arrival – Cell array of arrival structs (one per class) with fields: lambda, lambda0, lambda1, ma, Ia service: Cell array of service structs (one per class) with fields: mu, ST, St, tau_st, SerChoice K: Number of parallel queues fjInfo: Updated fjInfo with distribution info

Examples

[arrival, service, K, fjInfo] = fj_extract_params(sn, fjInfo)

References

fj_xmax_coxian(K, mu1, mu2, q)

Expected maximum of K i.i.d. two-stage Coxian variables

Exact expected maximum of K independent copies of X = T1 + B*T2, where T1 ~ Exp(mu1), T2 ~ Exp(mu2) and B is Bernoulli(q). The survival function is a two-term exponential mixture,

S(t) = A*exp(-mu1 t) + B*exp(-mu2 t), A = (1-q) + q*mu2/(mu2-mu1), B = -q*mu1/(mu2-mu1),

so expanding 1 - (1-S)^K binomially and integrating term by term gives

E[Y_K] = sum_{j=1..K} (-1)^(j+1) binom(K,j) sum_{i=0..j} binom(j,i) A^(j-i) B^i / ((j-i)*mu1 + i*mu2).

When the two stage rates coincide the mixture degenerates into S(t) = (1 + q*mu*t)*exp(-mu*t) and the same expansion is carried out with the polynomial integrals integral t^m exp(-j mu t) dt = m!/(j mu)^(m+1), which the implementation selects automatically.

Parameters:

K – Number of branches (positive integer)
mu1 – Rate of the first stage (positive)
mu2 – Rate of the second stage (positive)
q – Probability that the second stage is visited, in [0,1]

Returns:

Xmax – Exact expected maximum of the K branches m1: Mean of a single branch, 1/mu1 + q/mu2 c2: Squared coefficient of variation of a single branch

Examples

Xmax = fj_xmax_coxian(K, mu1, mu2, q)
[Xmax, m1, c2] = fj_xmax_coxian(K, mu1, mu2, q)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Section 4.5, Eqs. (38)-(40).

Original: P. M. Chen, D. Towsley, “A Performance Evaluation of RAID Architectures”, IEEE Trans. Computers 45(10), 1996.

mexify_fj: Create configuration object of class ‘coder.CodeConfig’.

fj_amva(D, P, M, Z)

Mean value analysis of a closed network of fork-join subnetworks

Population-by-population recursion for a single-class closed queueing network whose stations are P(n)-way fork-join subnetworks with per-visit demand D(n). The response time of a subnetwork inflates the arrival-instant queue length by the harmonic number of its fork degree, which is the Varki bound on the residence time of a parallel subsystem:

R_n(m) = D_n * [ H_{P_n} + Q_n(m-1) ] X(m) = m / (Z + sum_n R_n(m)) Q_n(m) = X(m) * R_n(m)

started from Q_n(0) = 0. Setting P(n) = 1 for every n recovers the exact single-class mean value analysis of Reiser and Lavenberg, because H_1 = 1. For P(n) > 1 the recursion is an approximation whose per-subnetwork residence time is an upper bound in the sense of Varki.

Parameters:

D – Vector of per-visit service demands, one per subnetwork
P – Vector of fork degrees, one per subnetwork (P(n) >= 1)
M – Number of circulating jobs (positive integer)
Z – Think time (optional, default 0)

Returns:

R – Vector of mean residence times, one per subnetwork Q: Vector of mean queue lengths, one per subnetwork X: System throughput U: Vector of utilizations of the busiest queue of each subnetwork

Examples

[R, Q, X, U] = fj_amva(D, P, M)
[R, Q, X, U] = fj_amva(D, P, M, Z)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014, Eqs. (68)-(69).

Original: E. Varki, “Mean Value Technique for Closed Fork-Join Networks”, ACM SIGMETRICS, 1999; G. Casale, R. R. Muntz, G. Serazzi, “Geometric Bounds”, IEEE Trans. Computers 57(6), 2008.

fj_xmax_hyperexp(K, p1, mu1, mu2)

Expected maximum of K i.i.d. Hyperexponential-2 random variables

Computes the expected value of the maximum of K independent and identically distributed Hyperexponential random variables with two branches.

The pdf of H2 is: f(t) = p1*mu1*exp(-mu1*t) + p2*mu2*exp(-mu2*t) where p1 + p2 = 1 and p1, p2 > 0.

The expected maximum is: X_K^max = sum_{n=1}^{K} (-1)^{n+1} * sum_{m=0}^{n} C(n,m) * p1^m * p2^{n-m} / (m*mu1 + (n-m)*mu2)

The H2 distribution has CV > 1, making it useful for modeling high-variability service times.

Parameters:

K – Number of parallel servers (positive integer)
p1 – Probability of branch 1 (0 < p1 < 1)
mu1 – Rate of branch 1
mu2 – Rate of branch 2

Returns:

Xmax – Expected maximum service time X_K^max

Examples

Xmax = fj_xmax_hyperexp(K, p1, mu1, mu2)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Page 17:13, Eq. (15).

sn_join_siblings(sn, joinIdx, r)

N = SN_JOIN_SIBLINGS(SN, JOINIDX, R)

Number of sibling tasks forked per parent job of class R on the fork-join pair that ends at the Join node JOINIDX. Siblings are counted at the FORK, as the simulation engines count them.

THE COUNT IS PER LINK, not the out-degree times a node-wide scalar. A fork carries a VARIABLE FORKING LEVEL: setTasksPerLink(class, n [, dest]) sets one link of one class, setTasksPerLinkDistribution makes the degree a draw, and setBranchProb makes a link taken with probability < 1. All three land in sn.nodeparam{f}.fanOutLink and .fanOutProb, both (nnodes x nclasses) and indexed by DESTINATION NODE, with the DISTRIBUTION case storing its mean. So

N = sum_d fanOutProb(d,r) * fanOutLink(d,r)

which is the EXPECTED sibling count, and reduces to out-degree times the node-wide fanOut on a model that sets none of the three (every connected link then carries fanOutLink = tasksPerLink and fanOutProb = 1).

Falls back to that older product when fanOutLink is absent – it is built by refreshLocalVars, which runs after refreshCapacity – and to the join’s in-degree when the matched fork cannot be identified.

N is what a quorum is measured against: the join fires on the k-th of N siblings and the remaining N-k are discarded when they arrive.

fj_tail_forktail(ET, VT, K, p, P)

ForkTail black-box tail-latency approximation for fork-join requests

Approximates the p-th percentile of the response time of a request that forks into K parallel tasks and joins on the last of them, from the mean and variance of the per-branch task response times alone. Each branch is treated as a black box: its task response time is fitted by a generalized exponential law

F_T(x) = (1 - exp(-x/beta))^alpha

whose two parameters are matched on the branch mean and variance,

E[T] = beta*(psi(alpha+1) - psi(1)), V[T] = beta^2*(psi’(1) - psi’(alpha+1)),

and the request response time is the maximum over the branches, taken as the product of the branch CDFs (exact only for independent branches):

F_X(x) = prod_i (1 - exp(-x/beta_i))^alpha_i, x_p = F_X^{-1}(p).

In the homogeneous case this inverts in closed form, x_p = -beta*log(1 - p^(1/(K*alpha))).

When the fanout itself is random, i.e. a request spawns K_i tasks with probability P_i (a service whose requests touch different numbers of shards), the request law is the mixture

F_X(x) = sum_i P_i * (1 - exp(-x/beta))^(K_i*alpha),

which is inverted numerically.

The approximation rests on the central limit theorem for G/G/m queues in heavy traffic, so it is a HIGH-LOAD result: the reference reports errors within 20% and 15% at 80% and 90% utilization respectively, and makes no claim at low load, where the tail is dominated by the service law rather than by queueing and the branch dependence is strongest. Use fj_is_homogeneous plus the FJ_codes route when the model is in the homogeneous MAP/PH/1 class, which is more accurate there; ForkTail covers the heterogeneous branches and mixed service laws that route rejects.

Parameters:

ET – Mean task response time: a scalar (homogeneous branches) or a vector with one entry per branch
VT – Variance of the task response time, same shape as ET
K – Number of branches, or a vector of distinct fanouts when the fanout is random; ignored when ET is a vector (default 1)
p – Percentile, either a fraction in (0,1) or a percentage in (0,100) (default 99)
P – Probabilities of the fanouts in K, required when K is a vector of more than one entry

Returns:

xp – Predicted p-th percentile of the request response time alpha: Fitted shape parameter(s) of the generalized exponential beta: Fitted scale parameter(s)

Examples

[xp, alpha, beta] = fj_tail_forktail(ET, VT, K, p)

References

M. Nguyen, S. Alesawi, N. Li, H. Che, H. Jiang, “ForkTail: A Black-Box Fork-Join Tail Latency Prediction Model for User-Facing Datacenter Workloads”, ACM HPDC 2018, pp. 206-217.

fj_branch_members(lqn, joinaidx)

FJ_BRANCH_MEMBERS Activities belonging to each branch of an AND-join.

MEMBERS = FJ_BRANCH_MEMBERS(LQN, JOINAIDX) returns a cell array with one entry per branch feeding the AND-join activity JOINAIDX. Each entry lists the global activity indices lying on that branch, from its head (the activity spawned by the AND-fork) to its tail (the immediate predecessor of the join).

A branch is recovered by walking backwards from each immediate predecessor of the join until an activity marked POST_AND is reached, that activity being the branch head spawned by the fork. Branches between a fork and its join are disjoint paths, so the walk is unambiguous.

fj_xmax_approx(K, mu_X, sigma_X, dist_type)

General approximation for expected maximum of K random variables

Approximates the expected value of the maximum of K i.i.d. random variables using the mean-variance approximation from David [1970].

X_K^max ≈ mu_X + sigma_X * G(K)

where G(K) depends on the distribution type: - Exponential: G(K) = H_K - 1 (where H_K is the K-th harmonic number) - Uniform: G(K) = sqrt(3) * (K-1) / (K+1) - EVD: G(K) = sqrt(6) * ln(K) / pi - General: G(K) <= (K-1) / sqrt(2K-1) (upper bound)

Parameters:

K – Number of random variables (positive integer)
mu_X – Mean of the distribution
sigma_X – Standard deviation of the distribution
dist_type – Distribution type: ‘exp’, ‘uniform’, ‘evd’, ‘bound’ (default: ‘exp’)

Returns:

Xmax – Approximate expected maximum X_K^max GK: The G(K) factor used

Examples

Xmax = fj_xmax_approx(K, mu_X, sigma_X, dist_type)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Eq. (22) and (23) on page 17:16.

Original: H.A. David, “Order Statistics”, Wiley, 1970.

sn_join_quorum(sn, joinIdx, r, nbranches)

K = SN_JOIN_QUORUM(SN, JOINIDX, R, NBRANCHES)

Number of sibling tasks the Join node JOINIDX waits for in class R, given that NBRANCHES of them are forked. A standard join, an absent declaration, a non-positive quorum and a quorum that is not smaller than the sibling count all return NBRANCHES, i.e. the ordinary AND-join: those are the four ways a join fires only when every sibling has arrived.

The count is the one the simulation engines apply (SolverLDES fixes it at FORK time and discards the stragglers when they reach the join), so an analytical solver reading it here charges the same synchronisation event.

fj_dist2fj(lineMAP, distType, sn, ist, r)

Convert LINE PH/MAP distribution to FJ_codes format

Converts LINE’s MAP representation {D0, D1} to FJ_codes format. For arrivals: returns (lambda, lambda0, lambda1, ma, Ia) For service: returns (mu, ST, St, tau_st, SerChoice)

Supported distributions: - Exponential (Exp) - 2-phase Hyperexponential (HyperExp with 2 phases) - 2-phase Erlang (Erlang with 2 phases) - 2-phase MAP (MAP with 2 states)

Parameters:

lineMAP – LINE MAP representation {D0, D1}
distType – ‘arrival’ or ‘service’
sn – Network structure (for getting rate and process type)
ist – Station index
r – Class index

Returns:

fjDist – FJ_codes compatible structure

Examples

fjDist = fj_dist2fj(lineMAP, distType, sn, ist, r)

References

fj_xmax_exp(K, mu)

Expected maximum of K i.i.d. exponential random variables

Computes the expected value of the maximum of K independent and identically distributed exponential random variables with rate mu (mean service time x = 1/mu).

X_K^max = H_K / mu = H_K * x

where H_K is the K-th Harmonic number.

This is a fundamental quantity in Fork-Join analysis. For a split-merge (SM) queueing system, the maximum throughput is lambda_K^SM = 1 / X_K^max.

For large K, X_K^max ≈ (ln(K) + gamma) / mu, where gamma ≈ 0.57721 is the Euler-Mascheroni constant.

Parameters:

K – Number of parallel servers (positive integer)
mu – Service rate (mean service time is 1/mu)

Returns:

Xmax – Expected maximum service time X_K^max = H_K / mu

Examples

Xmax = fj_xmax_exp(K, mu)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Page 17:2 and Eq. (10) on page 17:11.

sn_join_droprate(sn, TN, AN)

DROPRATEJOIN = SN_JOIN_DROPRATE(SN, TN, AN)

Rate at which sibling tasks are discarded at each Join, as an (nstations x nclasses) matrix that is zero away from the Join rows.

A Join is the one station where the loss identity ArvR - Tput does NOT hold, because the two rates are in different units: AN counts the SIBLINGS offered to the join (N per parent job) while TN counts the PARENT jobs released by it (one per synchronisation). Reading ArvR - Tput there reports (N-1)/N of the offered traffic as lost at every join, standard joins included, when a standard join loses nothing at all.

The siblings a join actually consumes are K per synchronisation, where K is the quorum (K = N on a standard join), so

DropRateJoin = max(0, AN - K*TN)

which is 0 for a standard join and (N-K)*TN for a quorum, the rate at which the stragglers of an already-fired parent are discarded on arrival.

This is the DERIVED value, exact given TN and AN. A solver that MEASURES the discards on its own sample path (SolverLDES) reports its own; see NetworkSolver.getAvgLossTable.

fj_tail_ordstat(ET, VT, K, p, kreq)

Tail latency of a k-of-n (quorum) fork-join request

Approximates the p-th percentile of the response time of a request that forks into N parallel tasks and joins on the K-th of them, from the mean and variance of the per-branch task response times alone. K = N is the ordinary AND-join and reproduces FJ_TAIL_FORKTAIL exactly; K = 1 is the first completion.

Each branch is treated as the same black box FJ_TAIL_FORKTAIL uses: its task response time is fitted by a generalized exponential law

F_i(x) = (1 - exp(-x/beta_i))^alpha_i

matched on the branch mean and variance by GE_FIT. The request completes once K of the N branches have, so its law is the K-th ORDER STATISTIC of independent, not identically distributed branch times,

F_X(x) = P(at least K of the N branches are done by x),

which is the upper tail of a Poisson-binomial with success probabilities F_i(x). It is evaluated by the standard convolution recurrence, which adds no cancellation, and inverted by bisection. With homogeneous branches the recurrence collapses to the regularized incomplete beta function I_{F(x)}(K, N-K+1), used directly.

At K = N both routes reduce TERM BY TERM to prod_i F_i(x), the product of the branch CDFs that FJ_TAIL_FORKTAIL inverts, so a full join evaluates exactly as it did before this function existed.

BRANCH INDEPENDENCE is assumed, as in ForkTail: the branches of one request are positively correlated through their shared arrival instant, so the true quorum percentile is somewhat larger than this one. The same heavy-traffic caveat applies, see FJ_TAIL_FORKTAIL.

Parameters:

ET – Mean task response time: a scalar (homogeneous branches) or a vector with one entry per branch
VT – Variance of the task response time, same shape as ET
K – Number of branches; used only when ET is a scalar (default 1)
p – Percentile, either a fraction in (0,1) or a percentage in (0,100) (default 99)
kreq – Quorum: the join fires on the kreq-th branch (default: every branch)

Returns:

xp – Predicted p-th percentile of the request response time alpha: Fitted shape parameter(s) of the generalized exponential beta: Fitted scale parameter(s)

Examples

[xp, alpha, beta] = fj_tail_ordstat(ET, VT, K, p, kreq)

References

M. Nguyen, S. Alesawi, N. Li, H. Che, H. Jiang, “ForkTail: A Black-Box Fork-Join Tail Latency Prediction Model for User-Facing Datacenter Workloads”, ACM HPDC 2018, pp. 206-217, for the branch law. A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys 47(2), Article 17, 2014, Sec. 3, for the quorum.

fj_sm_tput(K, mu)

Maximum throughput for Split-Merge queueing system

Computes the maximum throughput for a K-way Split-Merge (SM) queueing system with exponential service times.

In an SM system, all tasks of a request must complete before the next request can be issued. The maximum throughput is:

lambda_K^SM = 1 / X_K^max = mu / H_K

where X_K^max = H_K / mu is the expected maximum service time.

For comparison, the maximum throughput of a K-way F/J system is lambda_K^FJ = mu, which exceeds lambda_K^SM by a factor H_K.

Parameters:

K – Number of parallel servers (positive integer)
mu – Service rate at each server

Returns:

lambda_max – Maximum throughput lambda_K^SM = mu / H_K

Examples

lambda_max = fj_sm_tput(K, mu)

References

A. Thomasian, “Analysis of Fork/Join and Related Queueing Systems”, ACM Computing Surveys, Vol. 47, No. 2, Article 17, July 2014. Page 17:2.