@@ -7,11 +7,10 @@ import (
77 "os"
88 "strconv"
99
10- "github.com/dreading/gospecfunc/bessel"
1110 "github.com/filecoin-project/lotus/build"
11+ skellampmf "github.com/rvagg/go-skellam-pmf"
1212 "github.com/urfave/cli/v2"
1313 "golang.org/x/exp/constraints"
14- "gonum.org/v1/gonum/stat/distuv"
1514)
1615
1716var finalityCmd = & cli.Command {
@@ -25,6 +24,10 @@ var finalityCmd = &cli.Command{
2524 & cli.StringFlag {
2625 Name : "input" ,
2726 },
27+ & cli.IntFlag {
28+ Name : "target" ,
29+ Usage : "target epoch for which finality is calculated" ,
30+ },
2831 },
2932 ArgsUsage : "[inputFile]" ,
3033 Action : func (cctx * cli.Context ) error {
@@ -49,14 +52,15 @@ var finalityCmd = &cli.Command{
4952 return err
5053 }
5154
52- blocksPerEpoch := 5.0 // Expected number of blocks per epoch
53- byzantineFraction := 0.3 // Upper bound on the fraction of malicious nodes in the network
54- currentEpoch := len (chain ) - 1 // Current epoch (end of history)
55- targetEpoch := currentEpoch - 30 // Target epoch for which finality is calculated
55+ blocksPerEpoch := 5.0 // Expected number of blocks per epoch
56+ byzantineFraction := 0.3 // Upper bound on the fraction of malicious nodes in the network
57+ currentEpoch := len (chain ) - 1 // Current epoch (end of history)
58+ // targetEpoch := currentEpoch - 30 // Target epoch for which finality is calculated
59+ targetEpoch := cctx .Int ("target" )
5660
5761 finality := FinalityCalcValidator (chain , blocksPerEpoch , byzantineFraction , currentEpoch , targetEpoch )
5862
59- fmt .Fprintf (cctx .App .Writer , "Finality probability: %f \n " , finality )
63+ fmt .Fprintf (cctx .App .Writer , "Finality=%v @ %d for chain len=%d \n " , finality , targetEpoch , currentEpoch )
6064
6165 return nil
6266 },
@@ -88,15 +92,14 @@ func FinalityCalcValidator(chain []int, blocksPerEpoch float64, byzantineFractio
8892 sumExpectedAdversarialBlocksI += rateMaliciousBlocks
8993 sumChainBlocksI += chain [i - 1 ]
9094 // Poisson(k=k, lambda=sum(f*e))
91- prLi := distuv. Poisson { Lambda : sumExpectedAdversarialBlocksI }. Prob ( float64 (k + sumChainBlocksI ))
95+ prLi := poissonProb ( sumExpectedAdversarialBlocksI , float64 (k + sumChainBlocksI ))
9296 prL [k ] = math .Max (prL [k ], prLi )
9397
94- // Break if prL[k] becomes negligible
95- if k > 1 && prL [k ] < negligibleThreshold && prL [k ] < prL [k - 1 ] {
96- maxKL = k
97- prL = prL [:k + 1 ]
98- break
99- }
98+ }
99+ if k > 1 && prL [k ] < negligibleThreshold && prL [k ] < prL [k - 1 ] {
100+ maxKL = k
101+ prL = prL [:k + 1 ]
102+ break
100103 }
101104 }
102105
@@ -108,7 +111,7 @@ func FinalityCalcValidator(chain []int, blocksPerEpoch float64, byzantineFractio
108111
109112 // Calculate Pr(B=k) for each value of k
110113 for k := 0 ; k <= maxKB ; k ++ {
111- prB [k ] = distuv. Poisson { Lambda : float64 (currentEpoch - targetEpoch ) * rateMaliciousBlocks }. Prob ( float64 (k ))
114+ prB [k ] = poissonProb ( float64 (currentEpoch - targetEpoch )* rateMaliciousBlocks , float64 (k ))
112115
113116 // Break if prB[k] becomes negligible
114117 if k > 1 && prB [k ] < negligibleThreshold && prB [k ] < prB [k - 1 ] {
@@ -119,11 +122,11 @@ func FinalityCalcValidator(chain []int, blocksPerEpoch float64, byzantineFractio
119122 }
120123
121124 // Compute M
122- prHgt0 := 1 - distuv. Poisson { Lambda : rateHonestBlocks }. Prob ( 0 )
125+ prHgt0 := 1 - poissonProb ( rateHonestBlocks , 0 )
123126
124127 expZ := 0.0
125128 for k := 0 ; k < int (4 * blocksPerEpoch ); k ++ {
126- pmf := distuv. Poisson { Lambda : rateMaliciousBlocks }. Prob ( float64 (k ))
129+ pmf := poissonProb ( rateMaliciousBlocks , float64 (k ))
127130 expZ += ((rateHonestBlocks + float64 (k )) / math .Pow (2 , float64 (k ))) * pmf
128131 }
129132
@@ -132,7 +135,7 @@ func FinalityCalcValidator(chain []int, blocksPerEpoch float64, byzantineFractio
132135 prM := make ([]float64 , maxKM + 1 )
133136 for k := 0 ; k <= maxKM ; k ++ {
134137 for i := maxIM ; i > 0 ; i -- {
135- probMI := SkellamPMF (k , float64 (i )* rateMaliciousBlocks , float64 (i )* ratePublicChain )
138+ probMI := skellampmf . SkellamPMF (k , float64 (i )* rateMaliciousBlocks , float64 (i )* ratePublicChain )
136139
137140 // Break if probMI becomes negligible
138141 if probMI < negligibleThreshold && probMI < prM [k ] {
@@ -185,6 +188,17 @@ func FinalityCalcValidator(chain []int, blocksPerEpoch float64, byzantineFractio
185188
186189 return math .Min (prError , 1.0 )
187190}
191+ func poissonProb (lambda float64 , x float64 ) float64 {
192+ return math .Exp (poissonLogProb (lambda , x ))
193+ }
194+
195+ func poissonLogProb (lambda float64 , x float64 ) float64 {
196+ if x < 0 || math .Floor (x ) != x {
197+ return math .Inf (- 1 )
198+ }
199+ lg , _ := math .Lgamma (math .Floor (x ) + 1 )
200+ return x * math .Log (lambda ) - lambda - lg
201+ }
188202
189203func sum [T constraints.Integer | constraints.Float ](s []T ) T {
190204 var total T
@@ -210,46 +224,3 @@ func min(a, b int) int {
210224 }
211225 return b
212226}
213-
214- // SkellamPMF calculates the probability mass function (PMF) of a Skellam distribution.
215- //
216- // The Skellam distribution is the probability distribution of the difference
217- // of two independent Poisson random variables.
218- //
219- // Arguments:
220- // * k - The difference of two Poisson random variables.
221- // * mu1 - The expected value of the first Poisson distribution.
222- // * mu2 - The expected value of the second Poisson distribution.
223- //
224- // Returns:
225- // * A float64 representing the PMF of the Skellam distribution at k.
226- func SkellamPMF (k int , mu1 float64 , mu2 float64 ) float64 {
227- // Based on https://github.com/jsoares/rusty-skellam/blob/main/src/lib.rs
228-
229- // Return NaN if parameters outside range
230- if math .IsNaN (mu1 ) || mu1 <= 0 || math .IsNaN (mu2 ) || mu2 <= 0 {
231- return math .NaN ()
232- }
233-
234- // Parameterise and compute the Modified Bessel function of the first kind
235- nu := float64 (k )
236- z := complex (2.0 * math .Sqrt (mu1 * mu2 ), 0 )
237- besselResult := bessel .I (nu , z )
238-
239- // Compute the pmf
240- return math .Exp (- (mu1 + mu2 )) * math .Pow (mu1 / mu2 , nu / 2.0 ) * real (besselResult )
241- }
242-
243- /*
244- func main() {
245- seed := rand.NewSource(1)
246- random := rand.New(seed)
247- chain := make([]int, 1000)
248- for i := range chain {
249- chain[i] = random.Intn(5) + 1
250- }
251-
252- errorProbability := FinalityCalcValidator(chain, 5.0, 0.3, 1000, 900)
253- fmt.Printf("Error probability: %f\n", errorProbability)
254- }
255- */
0 commit comments