diff --git a/CHANGELOG.md b/CHANGELOG.md --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,11 @@ # Revision history for `deltaq` +## 1.2.1.0 - 2026-11-05 + +### Added + +* Add functions `onConditionOfSuccess` and `residualTime` to `DeltaQ.PiecewisePolynomial`. + ## 1.2.0.0 — 2025-11-21 ### Changed diff --git a/deltaq.cabal b/deltaq.cabal --- a/deltaq.cabal +++ b/deltaq.cabal @@ -5,7 +5,7 @@ -- PVP summary: +-+------- breaking API changes -- | | +----- non-breaking API additions -- | | | +--- code changes with no API change -version: 1.2.0.0 +version: 1.2.1.0 synopsis: Framework for ∆Q System Development description: ∆Q System Development is a paradigm for developing distributed systems @@ -29,16 +29,16 @@ license: BSD-3-Clause license-file: LICENSE copyright: - Predictable Network Solutions Ltd., 2003-2025 + Predictable Network Solutions Ltd., 2003-2026 PLWORKZ R&D, 2025 author: Neil Davies, Heinrich Apfelmus -maintainer: neil.davies@pnsol.ccom +maintainer: neil.davies@pnsol.com extra-doc-files: CHANGELOG.md README.md images/hops.svg -tested-with: GHC ==8.10.7 || ==9.6.7 || ==9.10.3 +tested-with: GHC ==8.10.7 || ==9.6.7 || ==9.10.3 || ==9.12.4 common warnings ghc-options: -Wall @@ -55,10 +55,10 @@ build-depends: , base >=4.14.3.0 && <5 , Chart >=1.8 && <2 - , deepseq >=1.4.4.0 && <1.6 + , deepseq >=1.4.4.0 && <1.7 , lattices >=2.2 && <2.3 - , probability-polynomial >=1.0.1 && <1.1 - , diagrams-lib >=1.4.5 && <1.6 + , probability-polynomial >=1.0.2 && <1.2 + , diagrams-lib >=1.4.5 && <1.7 , diagrams-svg >=1.5 && <1.6 exposed-modules: @@ -82,7 +82,7 @@ , deltaq , hspec >=2.11.0 && <2.12 , probability-polynomial - , QuickCheck >=2.14 && <2.18 + , QuickCheck >=2.14 && <2.19 main-is: Spec.hs other-modules: diff --git a/src/DeltaQ/PiecewisePolynomial.hs b/src/DeltaQ/PiecewisePolynomial.hs --- a/src/DeltaQ/PiecewisePolynomial.hs +++ b/src/DeltaQ/PiecewisePolynomial.hs @@ -24,7 +24,9 @@ , meetsQTA , Moments (..) , moments + , onConditionOfSuccess , timeout + , residualTime -- * Internal , complexity @@ -57,7 +59,9 @@ import Numeric.Probability.Moments ( Moments (..) ) - +import Numeric.Function.Piecewise.Polynomial + ( integral + ) import qualified Data.Function.Class as Function import qualified Numeric.Function.Piecewise as Piecewise import qualified Numeric.Measure.Finite.Mixed as Measure @@ -276,6 +280,30 @@ {----------------------------------------------------------------------------- Operations + Conditional on success +------------------------------------------------------------------------------} +-- | Compute the probability (distribution) on the condition +-- that the outcome completes in a finite amount of time. +-- +-- Put differently: Rescale the cumulative distribution function +-- such that the probability of failure becomes zero. +-- +-- Return 'Nothing' if the outcome never completes, 'Just' otherwise. +-- +-- > onConditionOfSuccess never = Nothing +-- > ∀x. failure x < 1 ==> failure (onConditionOfSuccess x) = Just 0 +-- > ∀x. failure x < 1 ==> +-- > (snd . moments <$> onConditionOfSuccess x) = Just (snd (moments x)) +onConditionOfSuccess :: DQ -> Maybe DQ +onConditionOfSuccess (DQ m) + | success == 0 = Nothing + | otherwise = Just $ DQ conditional + where + success = Measure.total m + conditional = Measure.scale (1/success) m + +{----------------------------------------------------------------------------- + Operations Timeout ------------------------------------------------------------------------------} -- | Given a duration @dt@, @@ -312,3 +340,50 @@ afterCond | p == 1 = never | otherwise = DQ $ Measure.scale (1/(1-p)) after + +{----------------------------------------------------------------------------- + Operations + Residual Time +------------------------------------------------------------------------------} +{- | + For a renewal process, + `residualTime` computes the time needed to wait until the next renewal, + given that we have already waited an indefinite/random amount of time. + + The input argument is the probability distribution of completion times + of the renewal process. + The result is the residual (equilibrium) distribution of waiting times. + + The residual time is computed as + $G(t) = (1/mu) * integral from 0 to t of (1 - F(s)) ds$. + See also + . + + Note that the input distribution must be a valid probability distribution + (i.e. total probability = 1), otherwise the result is `Nothing`. + + If it has zero mean, then it must be a Dirac distribution at 0; considering limits, + the residual time is then also a Dirac distribution at 0. +-} +residualTime :: DQ -> Maybe DQ +residualTime f + | success /= 1 = Nothing + | mu == 0 = Just f + | otherwise = unsafeFromPositiveMeasure <$> Measure.fromDistribution g + where + (success, moments') = moments f + mu = mean moments' + df = distribution f -- unpack the DQ to the piecewise polynomial representation + oneMinusF = + Piecewise.fromAscPieces + $ oneFromZero + $ Piecewise.toAscPieces + $ Piecewise.mapPieces (Poly.constant 1 -) df + g = Piecewise.mapPieces (Poly.scale (1/mu)) $ integral oneMinusF + +-- a Piecewise is zero before its first piece, but 1 - F is 1 on [0, x1), +-- so add that piece explicitly if necessary +oneFromZero :: (Num a, Ord a) => [(a, Poly a)] -> [(a, Poly a)] +oneFromZero ps@((x, _) : _) + | x > 0 = (0, Poly.constant 1) : ps +oneFromZero ps = ps diff --git a/test/DeltaQ/PiecewisePolynomialSpec.hs b/test/DeltaQ/PiecewisePolynomialSpec.hs --- a/test/DeltaQ/PiecewisePolynomialSpec.hs +++ b/test/DeltaQ/PiecewisePolynomialSpec.hs @@ -34,6 +34,8 @@ , fromPositiveMeasure , meetsQTA , moments + , onConditionOfSuccess + , residualTime , timeout ) import Numeric.Probability.Moments @@ -56,6 +58,7 @@ , arbitrary , choose , chooseInteger + , counterexample , frequency , getSize , mapSize @@ -466,6 +469,21 @@ let ms = Moments{mean = t, variance = 0, skewness = 0, kurtosis = 1} in moments (wait t) === (1, ms) + describe "onConditionOfSuccess" $ do + it "never" $ withMaxSuccess 1 $ property $ + onConditionOfSuccess never === Nothing + + it "failure" $ property $ + \x -> + failure x < 1 + ==> (failure <$> onConditionOfSuccess x) === Just 0 + + it "moments" $ property $ + \x -> + failure x < 1 + ==> (snd . moments <$> onConditionOfSuccess x) + === Just (snd (moments x)) + describe "timeout" $ do it "never" $ property $ \dt -> @@ -491,6 +509,57 @@ \x dt -> let (within, p, after) = timeout dt x in choice p within (wait dt .>>. after) === x + + describe "residualTime" $ do + it "never" $ withMaxSuccess 1 $ property $ + residualTime never === Nothing + + it "failure" $ property $ + \(Probability p) x -> + p < 1 ==> residualTime (choice p x never) === Nothing + + it "success" $ property $ + \x -> + failure x == 0 + ==> (failure <$> residualTime x) === Just 0 + + it "wait 0" $ withMaxSuccess 1 $ property $ + residualTime (wait 0) === Just (wait 0) + + it "wait" $ property $ + \(Positive t) -> + residualTime (wait t) === Just (uniform 0 t) + + -- For a renewal process with inter-event time X and mean mu = E[X], + -- the residual time R satisfies E[R^k] = E[X^(k+1)] / ((k+1) * mu). + -- + -- E[X^2] follows exactly from the mean and variance, + -- but E[X^3] requires the standard deviation, which 'moments' + -- only approximates, so the second moment is compared + -- up to a relative tolerance. + it "first two moments" $ property $ + \x -> + let (_, Moments{mean = mu, variance = var, skewness = skew}) = + moments x + m2 = var + mu * mu + sigma = sqrt (fromRational var) :: Double + m3 = fromRational skew * sigma ^ (3 :: Int) + + fromRational (3 * mu * var + mu * mu * mu) + in (failure x == 0 && mu > 0) ==> + case residualTime x of + Nothing -> + counterexample "residualTime returned Nothing" False + Just r -> + let (_, Moments{mean = muR, variance = varR}) = + moments r + expected = m3 / (3 * fromRational mu) + actual = fromRational (varR + muR * muR) + in muR === m2 / (2 * mu) + .&&. counterexample + ("E[R^2]: expected " <> show expected + <> ", got " <> show actual) + (abs (actual - expected) + <= 1e-3 * abs expected) describe "complexity" $ do it "grows exponentially with .>>." $ withMaxSuccess 1 $ property $