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
99
100
101
102
103
104
105
106
107
108
109
|
"""
Miscellaneous utilities.
"""
__all__ = ["draw", "console"]
import numpy as np
import numpy.random as rng
try:
from numba import njit
except ImportError:
def njit(*args, **kw):
return lambda f: f
@njit(cache=True)
def draw(k, n):
"""
Select k things from a pool of n without replacement.
"""
# At k == n/4, an extra 0.15*k draws are needed to get k unique draws
# TODO: silently returns too few values if k > n
if k > n / 4:
result = rng.permutation(n)[:k]
else:
result = np.empty(k, np.int64)
for i in range(k):
# select an item not already selected
while True:
p = rng.randint(n)
for j in range(i):
# if there is a match then break out of the for-loop and
# go to the next while iteration, generating a new proposal
if j == p:
break
# unusual syntax: if you make it to the end of the for loop
# then there were no matches, so break out of the while loop
else:
break
result[i] = p
return result
def _check_uniform_draw():
"""
Draws from history should
"""
import pylab
k, n = 50, 400
counts = np.zeros(n * k)
idx = np.arange(k)
for _ in range(100000):
t = draw(k, n)
counts[k * t + idx] += 1
pylab.subplot(211)
pylab.pcolormesh(np.reshape(counts, (n, k)))
pylab.colorbar()
pylab.title("drawn number vs. draw position")
pylab.subplot(212)
pylab.hist(counts)
pylab.title("number of draws per (number,position) bin")
pylab.show()
def console():
"""
Start the python console with the local variables available.
console() should be the last thing in the file, after sampling and
showing the default plots.
"""
import os
import sys
# Hack for eclipse console: can't actually run ipython in the eclipse
# console and get it to plot, so instead guess whether we are in a
# console by checking if we are attached to a proper tty through stdin.
# For eclipse, just show the plots.
try:
tty = os.isatty(sys.stdin.fileno())
except Exception:
tty = False
if tty:
# Display outstanding plots and turn interactive on
from matplotlib import interactive
from matplotlib._pylab_helpers import Gcf
for fig in Gcf.get_all_fig_managers():
try: # CRUFT
fig.show()
except AttributeError:
fig.frame.Show()
interactive(True)
# Start an ipython shell with the caller's local variables
import IPython
symbols = sys._getframe(1).f_locals
ip = IPython.Shell.IPShell(user_ns=symbols)
ip.mainloop()
else:
# Not a tty; try doing show() anyway
import pylab
pylab.show()
|
