forked from nytlabs/streamtools
/
poisson.go
98 lines (91 loc) · 2.18 KB
/
poisson.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
package library
import (
"errors"
"github.com/nytlabs/streamtools/st/blocks" // blocks
"github.com/nytlabs/streamtools/st/util" // util
"math"
"math/rand"
)
// specify those channels we're going to use to communicate with streamtools
type Poisson struct {
blocks.Block
queryrule chan chan interface{}
inrule chan interface{}
inpoll chan interface{}
out chan interface{}
quit chan interface{}
}
// we need to build a simple factory so that streamtools can make new blocks of this kind
func NewPoisson() blocks.BlockInterface {
return &Poisson{}
}
// Setup is called once before running the block. We build up the channels and specify what kind of block this is.
func (b *Poisson) Setup() {
b.Kind = "Poisson"
b.Desc = "draws a random number from a Poisson distribution when polled"
b.inrule = b.InRoute("rule")
b.queryrule = b.QueryRoute("rule")
b.inpoll = b.InRoute("poll")
b.quit = b.Quit()
b.out = b.Broadcast()
}
// algorithm due to Knuth http://en.wikipedia.org/wiki/Poisson_distribution
/*
algorithm poisson random number (Knuth):
init:
Let L ← e−^λ, k ← 0 and p ← 1.
do:
k ← k + 1.
Generate uniform random number u in [0,1] and let p ← p × u.
while p > L.
return k − 1.
*/
func NewPoissonSampler(λ float64) func() int {
L := math.Exp(-λ)
r := rand.New(rand.NewSource(12345))
return func() int {
k := 0
p := 1.0
for {
k = k + 1
u := r.Float64()
p = p * u
if p <= L {
return k - 1
}
}
}
}
func (b *Poisson) Run() {
var err error
λ := 1.0
sampler := NewPoissonSampler(λ)
for {
select {
case ruleI := <-b.inrule:
// set a parameter of the block
rule, ok := ruleI.(map[string]interface{})
if !ok {
b.Error(errors.New("couldn't assert rule to map"))
}
λ, err = util.ParseFloat(rule, "Rate")
if err != nil {
b.Error(err)
}
sampler = NewPoissonSampler(λ)
case <-b.quit:
// quit the block
return
case <-b.inpoll:
// deal with a poll request
b.out <- map[string]interface{}{
"sample": float64(sampler()),
}
case c := <-b.queryrule:
// deal with a query request
c <- map[string]interface{}{
"Rate": λ,
}
}
}
}