@@ -7,11 +7,11 @@ import (
77 "os"
88 "strconv"
99
10- "github.com/dreading/gospecfunc/bessel"
11- "github.com/filecoin-project/lotus/build"
10+ skellampmf "github.com/rvagg/go-skellam-pmf"
1211 "github.com/urfave/cli/v2"
1312 "golang.org/x/exp/constraints"
14- "gonum.org/v1/gonum/stat/distuv"
13+
14+ "github.com/filecoin-project/lotus/build"
1515)
1616
1717var finalityCmd = & cli.Command {
@@ -25,6 +25,10 @@ var finalityCmd = &cli.Command{
2525 & cli.StringFlag {
2626 Name : "input" ,
2727 },
28+ & cli.IntFlag {
29+ Name : "target" ,
30+ Usage : "target epoch for which finality is calculated" ,
31+ },
2832 },
2933 ArgsUsage : "[inputFile]" ,
3034 Action : func (cctx * cli.Context ) error {
@@ -33,7 +37,7 @@ var finalityCmd = &cli.Command{
3337 if err != nil {
3438 return err
3539 }
36- defer file .Close ()
40+ defer func () { _ = file .Close () } ()
3741
3842 var chain []int
3943 scanner := bufio .NewScanner (file )
@@ -49,14 +53,15 @@ var finalityCmd = &cli.Command{
4953 return err
5054 }
5155
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
56+ blocksPerEpoch := 5.0 // Expected number of blocks per epoch
57+ byzantineFraction := 0.3 // Upper bound on the fraction of malicious nodes in the network
58+ currentEpoch := len (chain ) - 1 // Current epoch (end of history)
59+ // targetEpoch := currentEpoch - 30 // Target epoch for which finality is calculated
60+ targetEpoch := cctx .Int ("target" )
5661
5762 finality := FinalityCalcValidator (chain , blocksPerEpoch , byzantineFraction , currentEpoch , targetEpoch )
5863
59- fmt .Fprintf (cctx .App .Writer , "Finality probability: %f \n " , finality )
64+ _ , _ = fmt .Fprintf (cctx .App .Writer , "Finality=%v @ %d for chain len=%d \n " , finality , targetEpoch , currentEpoch )
6065
6166 return nil
6267 },
@@ -69,10 +74,10 @@ func FinalityCalcValidator(chain []int, blocksPerEpoch float64, byzantineFractio
6974 // Threshold at which the probability of an event is considered negligible
7075 const negligibleThreshold = 1e-25
7176
72- maxKL := 400 // Max k for which to calculate Pr(L=k)
73- maxKB := int (( currentEpoch - targetEpoch ) * int (blocksPerEpoch ) ) // Max k for which to calculate Pr(B=k)
74- maxKM := 400 // Max k for which to calculate Pr(M=k)
75- maxIM := 100 // Maximum number of epochs for the calculation (after which the pr become negligible)
77+ maxKL := 400 // Max k for which to calculate Pr(L=k)
78+ maxKB := ( currentEpoch - targetEpoch ) * int (blocksPerEpoch ) // Max k for which to calculate Pr(B=k)
79+ maxKM := 400 // Max k for which to calculate Pr(M=k)
80+ maxIM := 100 // Maximum number of epochs for the calculation (after which the pr become negligible)
7681
7782 rateMaliciousBlocks := blocksPerEpoch * byzantineFraction // upper bound
7883 rateHonestBlocks := blocksPerEpoch - rateMaliciousBlocks // lower bound
@@ -88,15 +93,14 @@ func FinalityCalcValidator(chain []int, blocksPerEpoch float64, byzantineFractio
8893 sumExpectedAdversarialBlocksI += rateMaliciousBlocks
8994 sumChainBlocksI += chain [i - 1 ]
9095 // Poisson(k=k, lambda=sum(f*e))
91- prLi := distuv. Poisson { Lambda : sumExpectedAdversarialBlocksI }. Prob ( float64 (k + sumChainBlocksI ))
96+ prLi := poissonProb ( sumExpectedAdversarialBlocksI , float64 (k + sumChainBlocksI ))
9297 prL [k ] = math .Max (prL [k ], prLi )
9398
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- }
99+ }
100+ if k > 1 && prL [k ] < negligibleThreshold && prL [k ] < prL [k - 1 ] {
101+ maxKL = k
102+ prL = prL [:k + 1 ]
103+ break
100104 }
101105 }
102106
@@ -108,7 +112,7 @@ func FinalityCalcValidator(chain []int, blocksPerEpoch float64, byzantineFractio
108112
109113 // Calculate Pr(B=k) for each value of k
110114 for k := 0 ; k <= maxKB ; k ++ {
111- prB [k ] = distuv. Poisson { Lambda : float64 (currentEpoch - targetEpoch ) * rateMaliciousBlocks }. Prob ( float64 (k ))
115+ prB [k ] = poissonProb ( float64 (currentEpoch - targetEpoch )* rateMaliciousBlocks , float64 (k ))
112116
113117 // Break if prB[k] becomes negligible
114118 if k > 1 && prB [k ] < negligibleThreshold && prB [k ] < prB [k - 1 ] {
@@ -119,11 +123,11 @@ func FinalityCalcValidator(chain []int, blocksPerEpoch float64, byzantineFractio
119123 }
120124
121125 // Compute M
122- prHgt0 := 1 - distuv. Poisson { Lambda : rateHonestBlocks }. Prob ( 0 )
126+ prHgt0 := 1 - poissonProb ( rateHonestBlocks , 0 )
123127
124128 expZ := 0.0
125129 for k := 0 ; k < int (4 * blocksPerEpoch ); k ++ {
126- pmf := distuv. Poisson { Lambda : rateMaliciousBlocks }. Prob ( float64 (k ))
130+ pmf := poissonProb ( rateMaliciousBlocks , float64 (k ))
127131 expZ += ((rateHonestBlocks + float64 (k )) / math .Pow (2 , float64 (k ))) * pmf
128132 }
129133
@@ -132,7 +136,7 @@ func FinalityCalcValidator(chain []int, blocksPerEpoch float64, byzantineFractio
132136 prM := make ([]float64 , maxKM + 1 )
133137 for k := 0 ; k <= maxKM ; k ++ {
134138 for i := maxIM ; i > 0 ; i -- {
135- probMI := SkellamPMF (k , float64 (i )* rateMaliciousBlocks , float64 (i )* ratePublicChain )
139+ probMI := skellampmf . SkellamPMF (k , float64 (i )* rateMaliciousBlocks , float64 (i )* ratePublicChain )
136140
137141 // Break if probMI becomes negligible
138142 if probMI < negligibleThreshold && probMI < prM [k ] {
@@ -186,6 +190,18 @@ func FinalityCalcValidator(chain []int, blocksPerEpoch float64, byzantineFractio
186190 return math .Min (prError , 1.0 )
187191}
188192
193+ func poissonProb (lambda float64 , x float64 ) float64 {
194+ return math .Exp (poissonLogProb (lambda , x ))
195+ }
196+
197+ func poissonLogProb (lambda float64 , x float64 ) float64 {
198+ if x < 0 || math .Floor (x ) != x {
199+ return math .Inf (- 1 )
200+ }
201+ lg , _ := math .Lgamma (math .Floor (x ) + 1 )
202+ return x * math .Log (lambda ) - lambda - lg
203+ }
204+
189205func sum [T constraints.Integer | constraints.Float ](s []T ) T {
190206 var total T
191207 for _ , v := range s {
@@ -210,46 +226,3 @@ func min(a, b int) int {
210226 }
211227 return b
212228}
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