Package jline.solvers.fluid.moments

Class FluidClosures

java.lang.Object

jline.solvers.fluid.moments.FluidClosures

public final class FluidClosures
extends Object

Moment closures of the fluid drift, backing options.method='minnormal' of SolverFluid. Java twin of the MATLAB fluid_lld_scaling, fluid_min_closure, fluid_capacity_closure, fluid_share_closure and fluid_gps_share.

The mean-field ODEs of SolverFluid close the moment hierarchy at first order, replacing E[min(X,c)] by min(E[X],c) and the capacity share w_j X_j / sum_m w_m X_m by that ratio at the mean. Both are exact only where the map is locally linear, so their error peaks at the kink, i.e. exactly where the two means meet (rho ~ 1 at a queueing station). The functions here return instead the expectation under a normal marginal whose variance is produced by the covariance (Lyapunov) equation, so mean and covariance are solved self-consistently by SolverFluidMoments.

Every closure collapses to the first-order expression when its variance argument is zero, so callers share a single code path and the untouched methods are bit-identical.

See Also:

• FluidMomentTerms
• MinNormalAnalyzer

Nested Class Summary

Nested Classes

Modifier and Type	Class	Description
static final class	FluidClosures.ShareResult	A vector of capacity shares together with its Jacobian, and the joint-closure covariance Cov(S_j,N) and its Jacobian.
static final class	FluidClosures.ValueDeriv	A closure value together with its first two derivatives with respect to the mean.

Method Summary

Modifier and Type	Method	Description
static FluidClosures.ValueDeriv	capacityClosure(double n, double c, double s2, double[] lldrow, boolean isInf)	Moment closure of the station capacity term psi(X) and of its derivative, mirroring MATLAB fluid_capacity_closure.
static double[]	expansionWeight(double ratio)	How much of the second-order correction the series admits at this point, and d tau / d ratio, mirroring MATLAB local_expansion_weight.
static FluidClosures.ShareResult	gpsShare(Matrix xk, Matrix wk, Matrix vk, boolean wantJac)	Expected capacity share of a GPS station under a normal marginal, mirroring MATLAB fluid_gps_share.
static FluidClosures.ValueDeriv	lldScaling(double[] lldrow, double n)	Limited load-dependent rate scaling at a CONTINUOUS population, mirroring MATLAB fluid_lld_scaling.
static FluidClosures.ValueDeriv	minClosure(double n, double c, double s2)	Specialisation of minClosure(double, double, double, double, double) to a deterministic second argument.
static FluidClosures.ValueDeriv	minClosure(double n, double c, double s2, double vc, double covNc)	Min-normal moment closure of E[min(X,Y)] for jointly normal X, Y, mirroring MATLAB fluid_min_closure.
static double	normcdf(double z)	Standard normal cumulative distribution function.
static double	normpdf(double z)	Standard normal probability density function.
static void	projectRate(double[] r, double[] xb, boolean capped, double tot)	Project a jointly closed per-coordinate service share onto the set it has to live in, mirroring MATLAB local_project_rate: r >= 0, r <= xb where that bound applies, and sum(r) = tot.
static FluidClosures.ShareResult	shareClosure(Matrix x, Matrix wv, Matrix cov, boolean wantJac)	The share closure without the joint-closure covariance.
static FluidClosures.ShareResult	shareClosure(Matrix x, Matrix wv, Matrix cov, boolean wantJac, boolean wantCov)	Second-order closure of the capacity share of a sharing discipline, mirroring MATLAB fluid_share_closure.

Methods inherited from class java.lang.Object

clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait

Method Details

normcdf

public static double normcdf(double z)

Standard normal cumulative distribution function.

normpdf

public static double normpdf(double z)

Standard normal probability density function.

lldScaling

public static FluidClosures.ValueDeriv lldScaling(double[] lldrow,
 double n)

Limited load-dependent rate scaling at a CONTINUOUS population, mirroring MATLAB fluid_lld_scaling.

sn.lldscaling(i,:) tabulates the station rate multiplier at integer populations 1..lldlimit, and the discrete solvers read it as lldscaling(i, min(n, lldlimit)). The fluid state is continuous, so the table is read here by linear interpolation between consecutive entries, clamped to the first entry below n=1 and to the last entry above the table end, matching the clamping the CTMC already applies.

Parameters:

lldrow - rate multipliers at populations 1..lldlimit; null or empty means no load dependence and returns a = 1, da = 0

n - population, may be non-integer or negative

Returns:

alpha(n) and d alpha / d n, the latter zero on the clamped tails

minClosure

public static FluidClosures.ValueDeriv minClosure(double n,
 double c,
 double s2,
 double vc,
 double covNc)

Min-normal moment closure of E[min(X,Y)] for jointly normal X, Y, mirroring MATLAB fluid_min_closure.

This is the closure of Guenther, Stefanek and Bradley (EPEW/UKPEW 2012, LNCS 7587:32-47, eq. 4), in their general two-population form

   E[min(X,Y)] = E[X]*Phi((E[Y]-E[X])/th) + E[Y]*Phi((E[X]-E[Y])/th)
                 - th*phi((E[Y]-E[X])/th)
   th = (Var[X] - 2*Cov[X,Y] + Var[Y])^(1/2)
 

and the derivative with respect to E[X] is Phi((E[Y]-E[X])/th), the probability that X is the smaller of the two. SolverFluid only ever needs the specialisation Y = c, the deterministic server count: vc = 0 and covNc = 0 give th = sqrt(s2) and the expression collapses to the truncated-normal form. The general arguments are kept so this IS the published closure rather than one instance of it, and so a future state-dependent capacity needs no new derivation.

With th = 0 the expressions collapse to min(n,c) and to the indicator 1{n < c}, recovering the first-order closure exactly. The derivative at the kink is taken as 0, the right derivative of min(), which is the convention already implied by the strict inequality test in the closing rate factors.

Parameters:

n - mean of the first argument

c - mean of the second argument (server count when deterministic)

s2 - variance of the first argument

vc - variance of the second argument (0 when deterministic)

covNc - covariance of the two arguments

Returns:

E[min(X,Y)] and P(X < Y) under the normal marginal

minClosure

public static FluidClosures.ValueDeriv minClosure(double n,
 double c,
 double s2)

Specialisation of minClosure(double, double, double, double, double) to a deterministic second argument.

capacityClosure

public static FluidClosures.ValueDeriv capacityClosure(double n,
 double c,
 double s2,
 double[] lldrow,
 boolean isInf)

Moment closure of the station capacity term psi(X) and of its derivative, mirroring MATLAB fluid_capacity_closure.

Every scheduling branch of the closing rate factors scales the coordinates of a station by psi(n_i)/n_i, where psi is how much work the station clears at population n_i:

   psi(n) = min(n,c) * alpha(n)   at a queueing station
   psi(n) = n        * alpha(n)   at an infinite server
 

and alpha is the limited load-dependent scaling (1 when the station has none). This returns E[psi(X)] for X ~ Normal(n, s2) together with d/dn E[psi(X)]. s2 = 0 gives psi(n) itself, so the first-order closure is the same code path.

Without load dependence the expectation is the closed form of minClosure(double, double, double, double, double). With a tabulated alpha, lldScaling(double[], double) makes alpha piecewise linear on the integer lattice, so psi is piecewise QUADRATIC with breakpoints at the integers and at c, and the expectation is integrated segment by segment against the normal density using the truncated moments M0, M1, M2. The derivative is E[psi'(X)] by differentiation under the integral sign, valid because psi is Lipschitz.

The integration must be exact, NOT a fixed quadrature rule. Gauss-Hermite with fixed nodes applied to a piecewise-linear integrand does not smooth its kinks, it relocates them: the resulting estimate of E[psi] is itself piecewise linear in n, so its second derivative is zero almost everywhere and the 1/N refinement silently returns a null correction. That is what the segment-wise closed form below avoids.

psi is extended by zero below n = 0, the only physically admissible continuation, since a station holding no jobs clears no work.

Parameters:

n - mean population at the station

c - number of servers (ignored when isInf)

s2 - population variance, 0 for the first-order closure

lldrow - load-dependent scaling, null or empty when absent

isInf - true at an infinite-server station

Returns:

E[psi(X)] and d/dn E[psi(X)]

expansionWeight

public static double[] expansionWeight(double ratio)

How much of the second-order correction the series admits at this point, and d tau / d ratio, mirroring MATLAB local_expansion_weight.

Every second-order term of the share closure is a term of the series for E[1/v], whose successive terms are in the ratio Var(v)/v^2, so the truncation is meaningful below 1 and the terms GROW above it. Nothing in the algebra notices: at a near-empty station the corrections come back larger than the quantity they correct, and the drift that follows is not integrable.

One on [0,1], zero from 4 up, and the C^1 smoothstep between. Both ends matter. The lower one has to be EXACTLY one on the whole convergent region, so every model already inside it is bit-identical; the upper one has to be reached with a vanishing derivative, because the drift is integrated and a kink in it is what collapses the step size. The thresholds are the series, not a tuning: at ratio 1 successive terms stop shrinking, and at ratio 4 the standard deviation of v is twice its mean, where a non-negative v has essentially no mass near the point being expanded about.

Parameters:

ratio - Var(v)/v^2 at the point being expanded about

Returns:

{tau, d tau / d ratio}

projectRate

public static void projectRate(double[] r,
 double[] xb,
 boolean capped,
 double tot)

Project a jointly closed per-coordinate service share onto the set it has to live in, mirroring MATLAB local_project_rate: r >= 0, r <= xb where that bound applies, and sum(r) = tot.

THE JOINT CLOSURE IS AN EXPANSION AND CAN LEAVE THAT SET. r = s*psi + psi'*Cov(S,N) adds a term that sums to ZERO over the coordinates, so it moves mass between them and its entries can push one past either bound; the first-order share x_j/n_i*psi cannot, being x_j scaled by psi/n_i <= 1. Either breach ends the same way, because the integrator holds every coordinate non-negative: r_j > x_j drains coordinate j faster than it holds, the state goes negative and the clamp INJECTS mass.

THE UPPER BOUND HOLDS ONLY WITHOUT LOAD DEPENDENCE, which is what CAPPED selects. r is an expected NUMBER in service so r_j <= x_j, but psi(n) = min(n,c)*alpha(n) folds the load-dependent scaling into the same variable, and with alpha > 1 the first-order share itself exceeds x_j.

Clip, then move the residual onto the coordinates that still have slack in proportion to it, so sum(r) = tot survives and the station still clears what its capacity closure says it clears. It is a NO-OP whenever the expansion stayed inside the set, which is why models already inside it are bit-identical.

Parameters:

r - the jointly closed rates, projected in place

xb - the coordinate means, the upper bound when capped

capped - true without load dependence, where r_j <= x_j holds

tot - what the capacity closure says the station clears

shareClosure

public static FluidClosures.ShareResult shareClosure(Matrix x,
 Matrix wv,
 Matrix cov,
 boolean wantJac)

The share closure without the joint-closure covariance.

shareClosure

public static FluidClosures.ShareResult shareClosure(Matrix x,
 Matrix wv,
 Matrix cov,
 boolean wantJac,
 boolean wantCov)

Second-order closure of the capacity share of a sharing discipline, mirroring MATLAB fluid_share_closure.

A DPS station gives coordinate j the fraction S_j = w_j X_j / sum_m w_m X_m of its capacity, and PS is the same expression with unit weights. The first-order fluid closure evaluates that ratio at the mean, which is not E[S_j]: the map is a ratio, so Jensen's inequality biases it towards the coordinates carrying the LARGER weight. Writing u_j = w_j X_j and v = sum_m u_m, the delta method gives

   E[u_j/v] = mu_j/v - Cov(u_j,v)/v^2 + mu_j*Var(v)/v^3 + O(sigma^3)
 

evaluated at the means. The correction is exactly capacity-conserving: sum_j Cov(u_j,v) = Cov(v,v) = Var(v), so the two correction terms cancel in the sum and sum_j S_j = 1 identically, as it must for a work-conserving discipline. That identity is the invariant to check on any change here.

The expansion is local and fails when v is small against its own standard deviation, where the exact expectation is a Cauchy-like integral with no finite mean. There a raw share can come out negative; it is clipped at zero and the survivors renormalised, which preserves the conservation identity.

Parameters:

x - coordinate means of one station block (n x 1)

wv - per-coordinate weight, constant within a class (n x 1)

cov - covariance of the same coordinates (n x n); null or zero selects the first-order (plug-in) share

wantJac - whether to build the (n x n) Jacobian ds_j/dx_m

Returns:

the expected shares, summing to one, and their Jacobian

gpsShare

public static FluidClosures.ShareResult gpsShare(Matrix xk,
 Matrix wk,
 Matrix vk,
 boolean wantJac)

Expected capacity share of a GPS station under a normal marginal, mirroring MATLAB fluid_gps_share.

GPS divides the server by WEIGHT among the classes that are BACKLOGGED, and then equally among that class's own jobs (State.afterEventStation, GPS branch: cir = min(nir,1), share = w_r/(w*cir')). The share therefore depends on the backlog INDICATOR vector, not on the populations, and that is what makes GPS unreachable for a first-order closure: with continuous x_k > 0 every class is always backlogged, the indicator is identically one, and the share collapses to the constant w_r/sum_j w_j regardless of load. That constant is the heavy-traffic limit and is wrong everywhere else, so for GPS the second moment is not a correction, it is the entire mechanism.

The closure is an EXACT enumeration rather than an expansion. The share is piecewise CONSTANT over the 2^K backlog patterns, so

   E[S_r] = sum_{A contains r} P(backlog set = A) * w_r / sum_{j in A} w_j
 

carries no truncation error once the pattern probabilities are given. Those come from the marginals, P(N_k >= 1) = Phi((x_k - 1/2)/sigma_k) with the continuity correction for an integer population, multiplied as if the backlogs were independent. That independence is the one approximation here and it is not innocuous: in a closed network the station coordinates are NEGATIVELY correlated through population conservation, so the exact treatment would need multivariate-normal orthant probabilities.

The empty pattern contributes zero share, so sum_r E[S_r] = 1 - P(all classes empty) rather than 1. That is deliberate and is how the idle server is represented: GPS is single-server, so the backlog indicator plays the role that min(n,c) plays at a PS station, and no separate capacity term is applied by the caller.

Unlike the DPS ratio closure the expansion is NOT perturbative in sigma: as sigma_k -> 0 the probability tends to a step at x_k = 1 and its derivative phi(.)/sigma_k diverges, so the Jacobian stiffens at low variance. vk = 0 falls back to the hard indicator with zero derivative, which is the correct mean-field starting point for the outer iteration.

Parameters:

xk - per-class populations at the station (K x 1)

wk - per-class GPS weights, normalised internally (K x 1)

vk - per-class population variances; 0 selects the hard indicator

wantJac - whether to build the (K x K) Jacobian ds_r/dx_m

Returns:

the expected shares, summing to 1 - P(station empty), and their Jacobian