Here are the examples of the python api numpy.random.randn taken from open source projects. By voting up you can indicate which examples are most useful and appropriate.
175 Examples
3
Example 51
def __init__(self, nvis, nhid, iscale=0.1,
activation_fn=numpy.tanh,
params=None):
self.nvis = nvis
self.nhid = nhid
self.activation_fn = activation_fn
if params is None:
self.W = iscale * numpy.random.randn(nvis, nhid)
self.bias_vis = numpy.zeros(nvis)
self.bias_hid = numpy.zeros(nhid)
else:
self.W = params[0]
self.bias_vis = params[1]
self.bias_hid = params[2]
3
Example 52
Project: imitation Source File: nn.py
def test_standardizer():
D = 10
s = Standardizer(D, eps=0)
x_N_D = np.random.randn(200, D)
s.update(x_N_D)
x2_N_D = np.random.randn(300, D)
s.update(x2_N_D)
allx = np.concatenate([x_N_D, x2_N_D], axis=0)
assert np.allclose(s._mean_1_D.get_value()[0,:], allx.mean(axis=0))
assert np.allclose(s.get_stdev(), allx.std(axis=0))
print 'ok'
3
Example 53
Project: seaborn Source File: test_linearmodels.py
def test_three_point_colors(self):
x, y = np.random.randn(2, 3)
ax = lm.regplot(x, y, color=(1, 0, 0))
color = ax.collections[0].get_facecolors()
npt.assert_almost_equal(color[0, :3],
(1, 0, 0))
3
Example 54
def initialize_weights(self):
""" Randomly generate biases and weights for hidden layers.
Weights have a Gaussian distribution with mean 0 and
standard deviation 1 over the square root of the number
of weights connecting to the same neuron """
self.biases = [np.random.randn(y, 1) for y in self.layers[1:]]
self.weights = [np.random.randn(y, x)/np.sqrt(x) for x, y in zip(self.layers[:-1], self.layers[1:])]
3
Example 55
def predict(self, x_velocity, y_velocity, std ):
"""Predict the position of the point in the next frame.
Move the particles based on how the real system is predicted to behave.
The position of the point at the next time step is predicted using the
estimated velocity along X and Y axis and adding Gaussian noise sampled
from a distribution with MEAN=0.0 and STD=std. It is a linear model.
@param x_velocity the velocity of the object along the X axis in terms of pixels/frame
@param y_velocity the velocity of the object along the Y axis in terms of pixels/frame
@param std the standard deviation of the gaussian distribution used to add noise
"""
#To predict the position of the point at the next step we take the
#previous position and we add the estimated speed and Gaussian noise
self.particles[:, 0] += x_velocity + (np.random.randn(len(self.particles)) * std) #predict the X coord
self.particles[:, 1] += y_velocity + (np.random.randn(len(self.particles)) * std) #predict the Y coord
3
Example 56
Project: opendr Source File: test_geometry.py
def test_rodrigues(self):
from geometry import Rodrigues
rt = np.random.randn(3)
rt2 = rt + np.random.rand(3)*1e-5
foo1 = Rodrigues(rt = rt)
foo2 = Rodrigues(rt = rt2)
empirical = (foo2.r - foo1.r).flatten()
predicted = foo1.dr_wrt(foo1.rt).dot(rt2-rt)
self.assertTrue(np.max(np.abs(empirical - predicted)) < 1e-10)
3
Example 57
def setup(self):
# generate fake data
num_examples = 70
num_feats = 33
num_classes = 9 # EXCLUDING other class
self.x = np.random.randn(num_examples, num_feats)
self.cl = np.random.randint(0, num_classes, size=(num_examples))
self.l = np.zeros((num_examples, num_classes))
self.l[range(num_examples), self.cl] = 1
3
Example 58
def plot():
import matplotlib.pyplot as plt
import numpy as np
# Make plot with vertical (default) colorbar
fig = plt.figure()
ax = fig.add_subplot(111)
ax.hist(10 + 2 * np.random.randn(1000), label='men')
ax.hist(12 + 3 * np.random.randn(1000), label='women', alpha=0.5)
ax.legend()
return fig
3
Example 59
Project: data-analysis Source File: utils.py
def generate_measurements(x_0, dx, num_measurements, noise, acceleration=0):
data = []
for i in xrange(num_measurements):
data.append(x_0 + dx * i + np.random.randn() * noise)
dx += acceleration
return data
3
Example 60
def test_batch_size(self):
ds = downhill.Dataset([np.random.randn(40, 2)], batch_size=10, rng=4)
assert len(ds._slices) == 4
assert_size(ds, 0, 10)
assert_size(ds, 1, 10)
assert_size(ds, 2, 10)
assert_size(ds, 3, 10)
ds = downhill.Dataset([np.random.randn(40, 2)], batch_size=11, rng=4)
assert len(ds._slices) == 4
assert_size(ds, 0, 11)
assert_size(ds, 1, 11)
assert_size(ds, 2, 7)
assert_size(ds, 3, 11)
3
Example 61
def test_create(self, lgn):
series = random.randn(5, 100)
viz = lgn.plot(data={"series": series}, type='line')
assert isinstance(viz, Visualization)
assert hasattr(viz, 'id')
3
Example 62
def _reset(self):
"""
Reinitializes basis functions, iteration number and loss history.
"""
self.basis = np.random.randn(self.filterDim, self.outputDim)
self.basis /= np.sqrt(np.sum(self.basis ** 2, axis=0))
self._iteration = 0
self.losses = {}
3
Example 63
Project: py-videocore Source File: test_packunpack.py
def test_unpack_regA_float():
F = np.random.randn(16)
X = np.zeros((7, 16), dtype='uint32')
X[0] = unpack('16L', F.astype('float32'))
X[1] = unpack('16H', F.astype('float16'))
X[2] = unpack('16H', F.astype('float16'))
X[2] <<= 16
X[3:7] = np.array([getrandbits(32) for i in range(4*16)]).reshape(4, 16)
Y = run_code(unpack_regA_float, X)
assert all(X[0] == Y[0])
assert np.allclose(F, unpack('16f', Y[1]), rtol=1e-3)
assert np.allclose(F, unpack('16f', Y[2]), rtol=1e-3)
assert np.allclose(((X[3]>> 0)&0xff)/255.0, unpack('16f', Y[3]), rtol=1e-7)
assert np.allclose(((X[4]>> 8)&0xff)/255.0, unpack('16f', Y[4]), rtol=1e-7)
assert np.allclose(((X[5]>>16)&0xff)/255.0, unpack('16f', Y[5]), rtol=1e-7)
assert np.allclose(((X[6]>>24)&0xff)/255.0, unpack('16f', Y[6]), rtol=1e-7)
3
Example 64
def norm_weight(nin, nout=None, scale=0.01, ortho=True):
if nout is None:
nout = nin
if nout == nin and ortho:
W = ortho_weight(nin)
else:
W = scale * numpy.random.randn(nin, nout)
return W.astype('float32')
3
Example 65
Project: big_O Source File: test_big_o.py
def test_big_o(self):
desired = [(lambda n: [i for i in xrange(n*100)], compl.Linear),
(lambda n: 1., compl.Constant),
(lambda n: [i+j for i in xrange(n) for j in xrange(n)],
compl.Quadratic),
(lambda n: sorted(np.random.randn(n*100)),
compl.Linearithmic)]
for func, class_ in desired:
res_class, fitted = big_o.big_o(func, datagen.n_,
min_n=100, max_n=1000, n_repeats=5)
self.assertEqual(class_, res_class.__class__)
2
Example 66
Project: pyFFTW Source File: test_pyfftw_scipy_interface.py
def test_scipy_overwrite(self):
new_style_scipy_fftn = False
try:
scipy_fftn = scipy.signal.signaltools.fftn
scipy_ifftn = scipy.signal.signaltools.ifftn
except AttributeError:
scipy_fftn = scipy.fftpack.fftn
scipy_ifftn = scipy.fftpack.ifftn
new_style_scipy_fftn = True
a = pyfftw.empty_aligned((128, 64), dtype='complex128', n=16)
b = pyfftw.empty_aligned((128, 64), dtype='complex128', n=16)
a[:] = (numpy.random.randn(*a.shape) +
1j*numpy.random.randn(*a.shape))
b[:] = (numpy.random.randn(*b.shape) +
1j*numpy.random.randn(*b.shape))
scipy_c = scipy.signal.fftconvolve(a, b)
if new_style_scipy_fftn:
scipy.fftpack.fftn = scipy_fftpack.fftn
scipy.fftpack.ifftn = scipy_fftpack.ifftn
else:
scipy.signal.signaltools.fftn = scipy_fftpack.fftn
scipy.signal.signaltools.ifftn = scipy_fftpack.ifftn
scipy_replaced_c = scipy.signal.fftconvolve(a, b)
self.assertTrue(numpy.allclose(scipy_c, scipy_replaced_c))
if new_style_scipy_fftn:
scipy.fftpack.fftn = scipy_fftn
scipy.fftpack.ifftn = scipy_ifftn
else:
scipy.signal.signaltools.fftn = scipy_fftn
scipy.signal.signaltools.ifftn = scipy_ifftn
2
Example 67
Project: allantools Source File: noise.py
def iterpink(depth=20):
"""Generate a sequence of samples of pink noise.
pink noise generator
from http://pydoc.net/Python/lmj.sound/0.1.1/lmj.sound.noise/
Based on the Voss-McCartney algorithm, discussion and code examples at
http://www.firstpr.com.au/dsp/pink-noise/
depth: Use this many samples of white noise to calculate the output. A
higher number is slower to run, but renders low frequencies with more
correct power spectra.
Generates a never-ending sequence of floating-point values. Any continuous
set of these samples will tend to have a 1/f power spectrum.
"""
values = numpy.random.randn(depth)
smooth = numpy.random.randn(depth)
source = numpy.random.randn(depth)
sumvals = values.sum()
i = 0
while True:
yield sumvals + smooth[i]
# advance the index by 1. if the index wraps, generate noise to use in
# the calculations, but do not update any of the pink noise values.
i += 1
if i == depth:
i = 0
smooth = numpy.random.randn(depth)
source = numpy.random.randn(depth)
continue
# count trailing zeros in i
c = 0
while not (i >> c) & 1:
c += 1
# replace value c with a new source element
sumvals += source[i] - values[c]
values[c] = source[i]
2
Example 68
Project: pyFFTW Source File: test_pyfftw_call.py
def test_call_with_auto_input_alignment(self):
'''Test the class call with a keyword input update.
'''
input_array = (numpy.random.randn(*self.input_array.shape)
+ 1j*numpy.random.randn(*self.input_array.shape))
output_array = self.fft(
input_array=byte_align(input_array.copy(), n=16)).copy()
# Offset by one from 16 byte aligned to guarantee it's not
# 16 byte aligned
a = input_array
a__ = empty_aligned(numpy.prod(a.shape)*a.itemsize+1, dtype='int8',
n=16)
a_ = a__[1:].view(dtype=a.dtype).reshape(*a.shape)
a_[:] = a
# Just confirm that a usual update will fail
self.assertRaisesRegex(ValueError, 'Invalid input alignment',
self.fft.update_arrays, *(a_, self.output_array))
self.fft(a_, self.output_array)
self.assertTrue(numpy.alltrue(output_array == self.output_array))
# now try with a single byte offset and SIMD off
ar, ai = numpy.float32(numpy.random.randn(2, 257))
a = ar[1:] + 1j*ai[1:]
b = a.copy()
a_size = len(a.ravel())*a.itemsize
update_array = numpy.frombuffer(
numpy.zeros(a_size + 1, dtype='int8')[1:].data,
dtype=a.dtype).reshape(a.shape)
fft = FFTW(a, b, flags=('FFTW_UNALIGNED',))
# Confirm that a usual update will fail (it's not on the
# byte boundary)
self.assertRaisesRegex(ValueError, 'Invalid input alignment',
fft.update_arrays, *(update_array, b))
fft(update_array, b)
0
Example 69
def test_Helmholtz(ST2):
M = 4*N
kx = 12
points, weights = ST2.points_and_weights(M)
fj = np.random.randn(M)
f_hat = np.zeros(M)
if not ST2.__class__.__name__ == "ChebyshevTransform":
f_hat = ST2.fst(fj, f_hat)
fj = ST2.ifst(f_hat, fj)
if ST2.__class__.__name__ == "ShenDirichletBasis":
A = ADDmat(np.arange(M).astype(np.float))
B = BDDmat(np.arange(M).astype(np.float), ST2.quad)
s = slice(0, M-2)
elif ST2.__class__.__name__ == "ShenNeumannBasis":
A = ANNmat(np.arange(M).astype(np.float))
B = BNNmat(np.arange(M).astype(np.float), ST2.quad)
s = slice(1, M-2)
f_hat = np.zeros(M)
f_hat = ST2.fastShenScalar(fj, f_hat)
u_hat = np.zeros(M)
u_hat[s] = la.spsolve(A.diags()+kx**2*B.diags(), f_hat[s])
u1 = np.zeros(M)
u1 = ST2.ifst(u_hat, u1)
c = A.matvec(u_hat)+kx**2*B.matvec(u_hat)
c2 = np.dot(A.diags().toarray(), u_hat[s]) + kx**2*np.dot(B.diags().toarray(), u_hat[s])
assert np.allclose(c, f_hat)
assert np.allclose(c[s], c2)
H = Helmholtz(M, kx, ST2.quad, ST2.__class__.__name__ == "ShenNeumannBasis")
u0_hat = np.zeros(M)
u0_hat = H(u0_hat, f_hat)
u0 = np.zeros(M)
u0 = ST2.ifst(u0_hat, u0)
assert np.linalg.norm(u0 - u1) < 1e-12
# Multidimensional
f_hat = (f_hat.repeat(16).reshape((M, 4, 4))+1j*f_hat.repeat(16).reshape((M, 4, 4)))
kx = np.zeros((4, 4))+12
H = Helmholtz(M, kx, ST2.quad, ST2.__class__.__name__ == "ShenNeumannBasis")
u0_hat = np.zeros((M, 4, 4), dtype=np.complex)
u0_hat = H(u0_hat, f_hat)
u0 = np.zeros((M, 4, 4), dtype=np.complex)
u0 = ST2.ifst(u0_hat, u0)
#from IPython import embed; embed()
assert np.linalg.norm(u0[:, 2, 2].real - u1)/(M*16) < 1e-12
assert np.linalg.norm(u0[:, 2, 2].imag - u1)/(M*16) < 1e-12
0
Example 70
def test_noisy_1d():
f = KalmanFilter (dim_x=2, dim_z=1)
f.x = np.array([[2.],
[0.]]) # initial state (location and velocity)
f.F = np.array([[1.,1.],
[0.,1.]]) # state transition matrix
f.H = np.array([[1.,0.]]) # Measurement function
f.P *= 1000. # covariance matrix
f.R = 5 # state uncertainty
f.Q = 0.0001 # process uncertainty
fsq = SquareRootKalmanFilter (dim_x=2, dim_z=1)
fsq.x = np.array([[2.],
[0.]]) # initial state (location and velocity)
fsq.F = np.array([[1.,1.],
[0.,1.]]) # state transition matrix
fsq.H = np.array([[1.,0.]]) # Measurement function
fsq.P = np.eye(2) * 1000. # covariance matrix
fsq.R = 5 # state uncertainty
fsq.Q = 0.0001 # process uncertainty
measurements = []
results = []
zs = []
for t in range (100):
# create measurement = t plus white noise
z = t + random.randn()*20
zs.append(z)
# perform kalman filtering
f.update(z)
f.predict()
fsq.update(z)
fsq.predict()
assert abs(f.x[0,0] - fsq.x[0,0]) < 1.e-12
assert abs(f.x[1,0] - fsq.x[1,0]) < 1.e-12
# save data
results.append (f.x[0,0])
measurements.append(z)
p = f.P - fsq.P
print(f.P)
print(fsq.P)
for i in range(f.P.shape[0]):
assert abs(f.P[i,i] - fsq.P[i,i]) < 0.01
# now do a batch run with the stored z values so we can test that
# it is working the same as the recursive implementation.
# give slightly different P so result is slightly different
f.x = np.array([[2.,0]]).T
f.P = np.eye(2)*100.
m,c,_,_ = f.batch_filter(zs,update_first=False)
# plot data
if DO_PLOT:
p1, = plt.plot(measurements,'r', alpha=0.5)
p2, = plt.plot (results,'b')
p4, = plt.plot(m[:,0], 'm')
p3, = plt.plot ([0,100],[0,100], 'g') # perfect result
plt.legend([p1,p2, p3, p4],
["noisy measurement", "KF output", "ideal", "batch"], 4)
plt.show()
0
Example 71
Project: tensorly Source File: tucker_regression.py
def fit(self, X, y):
"""Fits the model to the data (X, y)
Parameters
----------
X : ndarray of shape (n_samples, N1, ..., NS)
tensor data
y : array of shape (n_samples)
labels associated with each sample
Returns
-------
self
"""
rng = check_random_state(self.random_state)
# Initialise randomly the weights
G = rng.randn(*self.weight_ranks)
W = []
for i in range(1, X.ndim): # First dimension of X = number of samples
W.append(np.random.randn(X.shape[i], G.shape[i - 1]))
# Norm of the weight tensor at each iteration
norm_W = []
for iteration in range(self.n_iter_max):
# Optimise modes of W
for i in range(len(W)):
phi = partial_tensor_to_vec(
np.dot(partial_unfold(X, i),
np.dot(kronecker(W, skip_matrix=i),
unfold(G, i).T)))
# Regress phi on y: we could call a package here, e.g. scikit-learn
inv_term = np.dot(phi.T, phi) + self.reg_W * np.eye(phi.shape[1])
W_i = vec_to_tensor(solve(inv_term, phi.T.dot(y)),
(X.shape[i + 1], G.shape[i]))
W[i] = W_i
phi = partial_tensor_to_vec(X).dot(kronecker(W))
G = vec_to_tensor(solve(phi.T.dot(phi) + self.reg_W * np.eye(phi.shape[1]), phi.T.dot(y)), G.shape)
weight_tensor_ = tucker_to_tensor(G, W)
norm_W.append(norm(weight_tensor_, 2))
# Convergence check
if iteration > 1:
weight_evolution = abs(norm_W[-1] - norm_W[-2]) / norm_W[-1]
if (weight_evolution <= self.tol):
if self.verbose:
print('\nConverged in {} iterations'.format(iteration))
break
self.weight_tensor_ = weight_tensor_
self.tucker_weight_ = (G, W)
self.vec_W_ = tucker_to_vec(G, W)
self.n_iterations_ = iteration + 1
self.norm_W_ = norm_W
return self
0
Example 72
Project: filterpy Source File: test_ukf.py
def test_circle():
from filterpy.kalman import KalmanFilter
from math import radians
def hx(x):
radius = x[0]
angle = x[1]
x = cos(radians(angle)) * radius
y = sin(radians(angle)) * radius
return np.array([x, y])
def fx(x, dt):
return np.array([x[0], x[1]+x[2], x[2]])
std_noise = .1
sp = JulierSigmaPoints(n=3, kappa=0.)
f = UKF(dim_x=3, dim_z=2, dt=.01, hx=hx, fx=fx, points=sp)
f.x = np.array([50., 90., 0])
f.P *= 100
f.R = np.eye(2)*(std_noise**2)
f.Q = np.eye(3)*.001
f.Q[0,0]=0
f.Q[2,2]=0
kf = KalmanFilter(dim_x=6, dim_z=2)
kf.x = np.array([50., 0., 0, 0, .0, 0.])
F = np.array([[1., 1., .5, 0., 0., 0.],
[0., 1., 1., 0., 0., 0.],
[0., 0., 1., 0., 0., 0.],
[0., 0., 0., 1., 1., .5],
[0., 0., 0., 0., 1., 1.],
[0., 0., 0., 0., 0., 1.]])
kf.F = F
kf.P*= 100
kf.H = np.array([[1,0,0,0,0,0],
[0,0,0,1,0,0]])
kf.R = f.R
kf.Q[0:3, 0:3] = Q_discrete_white_noise(3, 1., .00001)
kf.Q[3:6, 3:6] = Q_discrete_white_noise(3, 1., .00001)
measurements = []
results = []
zs = []
kfxs = []
for t in range (0,12000):
a = t / 30 + 90
x = cos(radians(a)) * 50.+ randn() * std_noise
y = sin(radians(a)) * 50. + randn() * std_noise
# create measurement = t plus white noise
z = np.array([x,y])
zs.append(z)
f.predict()
f.update(z)
kf.predict()
kf.update(z)
# save data
results.append (hx(f.x))
kfxs.append(kf.x)
#print(f.x)
results = np.asarray(results)
zs = np.asarray(zs)
kfxs = np.asarray(kfxs)
print(results)
if DO_PLOT:
plt.plot(zs[:,0], zs[:,1], c='r', label='z')
plt.plot(results[:,0], results[:,1], c='k', label='UKF')
plt.plot(kfxs[:,0], kfxs[:,3], c='g', label='KF')
plt.legend(loc='best')
plt.axis('equal')
0
Example 73
Project: pybrain Source File: cmaes.py
def _learnStep(self):
# Generate and evaluate lambda offspring
arz = randn(self.numParameters, self.batchSize)
arx = tile(self.center.reshape(self.numParameters, 1), (1, self.batchSize))\
+ self.stepSize * dot(dot(self.B, self.D), arz)
arfitness = zeros(self.batchSize)
for k in range(self.batchSize):
arfitness[k] = self._oneEvaluation(arx[:, k])
# Sort by fitness and compute weighted mean into center
arfitness, arindex = sorti(arfitness) # minimization
arz = arz[:, arindex]
arx = arx[:, arindex]
arzsel = arz[:, range(self.mu)]
arxsel = arx[:, range(self.mu)]
arxmut = arxsel - tile(self.center.reshape(self.numParameters, 1), (1, self.mu))
zmean = dot(arzsel, self.weights)
self.center = dot(arxsel, self.weights)
if self.storeAllCenters:
self._allCenters.append(self.center)
# Cumulation: Update evolution paths
self.stepPath = (1 - self.cuemStep) * self.stepPath \
+ sqrt(self.cuemStep * (2 - self.cuemStep) * self.muEff) * dot(self.B, zmean) # Eq. (4)
hsig = norm(self.stepPath) / sqrt(1 - (1 - self.cuemStep) ** (2 * self.numEvaluations / float(self.batchSize))) / self.chiN \
< 1.4 + 2. / (self.numParameters + 1)
self.covPath = (1 - self.cuemCov) * self.covPath + hsig * \
sqrt(self.cuemCov * (2 - self.cuemCov) * self.muEff) * dot(dot(self.B, self.D), zmean) # Eq. (2)
# Adapt covariance matrix C
self.C = ((1 - self.covLearningRate) * self.C # regard old matrix % Eq. (3)
+ self.covLearningRate * (1 / self.muCov) * (outer(self.covPath, self.covPath) # plus rank one update
+ (1 - hsig) * self.cuemCov * (2 - self.cuemCov) * self.C)
+ self.covLearningRate * (1 - 1 / self.muCov) # plus rank mu update
* dot(dot(arxmut, diag(self.weights)), arxmut.T)
)
# Adapt step size self.stepSize
self.stepSize *= exp((self.cuemStep / self.dampings) * (norm(self.stepPath) / self.chiN - 1)) # Eq. (5)
# Update B and D from C
# This is O(n^3). When strategy internal CPU-time is critical, the
# next three lines should be executed only every (alpha/covLearningRate/N)-th
# iteration, where alpha is e.g. between 0.1 and 10
self.C = (self.C + self.C.T) / 2 # enforce symmetry
Ev, self.B = eig(self.C) # eigen decomposition, B==normalized eigenvectors
Ev = real(Ev) # enforce real value
self.D = diag(sqrt(Ev)) #diag(ravel(sqrt(Ev))) # D contains standard deviations now
self.B = real(self.B)
# convergence is reached
if arfitness[0] == arfitness[-1] or (abs(arfitness[0] - arfitness[-1]) /
(abs(arfitness[0]) + abs(arfitness[-1]))) <= self.stopPrecision:
if self.verbose:
print("Converged.")
self.maxLearningSteps = self.numLearningSteps
# or diverged, unfortunately
if min(Ev) > 1e5:
if self.verbose:
print("Diverged.")
self.maxLearningSteps = self.numLearningSteps
0
Example 74
Project: pymanopt Source File: test_autograd.py
def setUp(self):
np.seterr(all='raise')
def f(x):
return (np.exp(np.sum(x[0]**2)) + np.exp(np.sum(x[1]**2)) +
np.exp(np.sum(x[2]**2)))
self.cost = f
n1 = self.n1 = 3
n2 = self.n2 = 4
n3 = self.n3 = 5
n4 = self.n4 = 6
n5 = self.n5 = 7
n6 = self.n6 = 8
self.y = y = (rnd.randn(n1), rnd.randn(n2, n3), rnd.randn(n4, n5, n6))
self.a = a = (rnd.randn(n1), rnd.randn(n2, n3), rnd.randn(n4, n5, n6))
self.correct_cost = f(y)
# CALCULATE CORRECT GRAD
g1 = 2 * y[0] * np.exp(np.sum(y[0] ** 2))
g2 = 2 * y[1] * np.exp(np.sum(y[1] ** 2))
g3 = 2 * y[2] * np.exp(np.sum(y[2] ** 2))
self.correct_grad = (g1, g2, g3)
# CALCULATE CORRECT HESS
# 1. VECTOR
Ymat = np.matrix(y[0])
Amat = np.matrix(a[0])
diag = np.eye(n1)
H = np.exp(np.sum(y[0] ** 2)) * (4 * Ymat.T.dot(Ymat) + 2 * diag)
# Then 'left multiply' H by A
h1 = np.array(Amat.dot(H)).flatten()
# 2. MATRIX
# First form hessian tensor H (4th order)
Y1 = y[1].reshape(n2, n3, 1, 1)
Y2 = y[1].reshape(1, 1, n2, n3)
# Create an m x n x m x n array with diag[i,j,k,l] == 1 iff
# (i == k and j == l), this is a 'diagonal' tensor.
diag = np.eye(n2 * n3).reshape(n2, n3, n2, n3)
H = np.exp(np.sum(y[1] ** 2)) * (4 * Y1 * Y2 + 2 * diag)
# Then 'right multiply' H by A
Atensor = a[1].reshape(1, 1, n2, n3)
h2 = np.sum(H * Atensor, axis=(2, 3))
# 3. Tensor3
# First form hessian tensor H (6th order)
Y1 = y[2].reshape(n4, n5, n6, 1, 1, 1)
Y2 = y[2].reshape(1, 1, 1, n4, n5, n6)
# Create an n1 x n2 x n3 x n1 x n2 x n3 diagonal tensor
diag = np.eye(n4 * n5 * n6).reshape(n4, n5, n6, n4, n5, n6)
H = np.exp(np.sum(y[2] ** 2)) * (4 * Y1 * Y2 + 2 * diag)
# Then 'right multiply' H by A
Atensor = a[2].reshape(1, 1, 1, n4, n5, n6)
h3 = np.sum(H * Atensor, axis=(3, 4, 5))
self.correct_hess = (h1, h2, h3)
self.backend = AutogradBackend()
0
Example 75
Project: filterpy Source File: test_ukf.py
def two_radar():
# code is not complete - I was using to test RTS smoother. very similar
# to two_radary.py in book.
import numpy as np
import matplotlib.pyplot as plt
from numpy import array
from numpy.linalg import norm
from numpy.random import randn
from math import atan2, radians
from filterpy.common import Q_discrete_white_noise
class RadarStation(object):
def __init__(self, pos, range_std, bearing_std):
self.pos = asarray(pos)
self.range_std = range_std
self.bearing_std = bearing_std
def reading_of(self, ac_pos):
""" Returns range and bearing to aircraft as tuple. bearing is in
radians.
"""
diff = np.subtract(self.pos, ac_pos)
rng = norm(diff)
brg = atan2(diff[1], diff[0])
return rng, brg
def noisy_reading(self, ac_pos):
rng, brg = self.reading_of(ac_pos)
rng += randn() * self.range_std
brg += randn() * self.bearing_std
return rng, brg
class ACSim(object):
def __init__(self, pos, vel, vel_std):
self.pos = asarray(pos, dtype=float)
self.vel = asarray(vel, dtype=float)
self.vel_std = vel_std
def update(self):
vel = self.vel + (randn() * self.vel_std)
self.pos += vel
return self.pos
dt = 1.
def hx(x):
r1, b1 = hx.R1.reading_of((x[0], x[2]))
r2, b2 = hx.R2.reading_of((x[0], x[2]))
return array([r1, b1, r2, b2])
pass
def fx(x, dt):
x_est = x.copy()
x_est[0] += x[1]*dt
x_est[2] += x[3]*dt
return x_est
vx, vy = 0.1, 0.1
f = UnscentedKalmanFilter(dim_x=4, dim_z=4, dt=dt, hx=hx, fx=fx, kappa=0)
aircraft = ACSim ((100,100), (vx*dt,vy*dt), 0.00000002)
range_std = 0.001 # 1 meter
bearing_std = 1/1000 # 1mrad
R1 = RadarStation ((0,0), range_std, bearing_std)
R2 = RadarStation ((200,0), range_std, bearing_std)
hx.R1 = R1
hx.R2 = R2
f.x = array([100, vx, 100, vy])
f.R = np.diag([range_std**2, bearing_std**2, range_std**2, bearing_std**2])
q = Q_discrete_white_noise(2, var=0.0002, dt=dt)
#q = np.array([[0,0],[0,0.0002]])
f.Q[0:2, 0:2] = q
f.Q[2:4, 2:4] = q
f.P = np.diag([.1, 0.01, .1, 0.01])
track = []
zs = []
for i in range(int(300/dt)):
pos = aircraft.update()
r1, b1 = R1.noisy_reading(pos)
r2, b2 = R2.noisy_reading(pos)
z = np.array([r1, b1, r2, b2])
zs.append(z)
track.append(pos.copy())
zs = asarray(zs)
xs, Ps, Pxz, pM, pP = f.batch_filter(zs)
ms, _, _ = f.rts_smoother(xs, Ps)
track = asarray(track)
time = np.arange(0,len(xs)*dt, dt)
plt.figure()
plt.subplot(411)
plt.plot(time, track[:,0])
plt.plot(time, xs[:,0])
plt.legend(loc=4)
plt.xlabel('time (sec)')
plt.ylabel('x position (m)')
plt.tight_layout()
plt.subplot(412)
plt.plot(time, track[:,1])
plt.plot(time, xs[:,2])
plt.legend(loc=4)
plt.xlabel('time (sec)')
plt.ylabel('y position (m)')
plt.tight_layout()
plt.subplot(413)
plt.plot(time, xs[:,1])
plt.plot(time, ms[:,1])
plt.legend(loc=4)
plt.ylim([0, 0.2])
plt.xlabel('time (sec)')
plt.ylabel('x velocity (m/s)')
plt.tight_layout()
plt.subplot(414)
plt.plot(time, xs[:,3])
plt.plot(time, ms[:,3])
plt.ylabel('y velocity (m/s)')
plt.legend(loc=4)
plt.xlabel('time (sec)')
plt.tight_layout()
plt.show()
0
Example 76
Project: MJHMC Source File: distributions.py
def cached_init_X(self):
""" Sets self.Xinit to cached (serialized) initial states for continuous-time samplers, generated by burn in
*For use with continuous-time samplers only*
:returns: None
:rtype: none
"""
distr_name = type(self).__name__
distr_hash = hash(self)
file_name = '{}_{}.pickle'.format(distr_name, distr_hash)
file_prefix = '{}/initializations'.format(package_path())
if file_name in os.listdir(file_prefix):
with open('{}/{}'.format(file_prefix, file_name), 'rb') as cache_file:
mjhmc_endpt, _, _, control_endpt = pickle.load(cache_file)
if self.mjhmc:
self.Xinit = mjhmc_endpt[:, :self.nbatch]
else:
self.Xinit = control_endpt[:, :self.nbatch]
else:
from mjhmc.misc.gen_mj_init import MAX_N_PARTICLES, cache_initialization
# modify this object so it can be used by gen_mj_init
old_nbatch = self.nbatch
self.nbatch = self.max_n_particles or MAX_N_PARTICLES
self.generation_instance = True
# must rebuild now that nbatch is changed back
if self.backend == 'tensorflow':
self.build_graph()
# start with biased initializations
# changes self.nbatch
try:
self.gen_init_X()
except NotImplementedError:
# completely arbitrary choice
self.Xinit = np.random.randn(self.ndims, self.nbatch)
#generate and cache fair initialization
cache_initialization(self)
# reconstruct this object using fair initialization
self.nbatch = old_nbatch
self.generation_instance = False
# must rebuild now that nbatch is changed back
if self.backend == 'tensorflow':
self.build_graph()
self.cached_init_X()
0
Example 77
Project: spectralDNS Source File: MKM.py
def initialize(solver, context):
# Initialize with pertubation ala perturbU (https://github.com/wyldckat/perturbU) for openfoam
U = context.U
X = context.X
params = config.params
Y = where(X[0]<0, 1+X[0], 1-X[0])
utau = params.nu * params.Re_tau
#Um = 46.9091*utau
Um = 56.*utau
Xplus = Y*params.Re_tau
Yplus = X[1]*params.Re_tau
Zplus = X[2]*params.Re_tau
duplus = Um*0.2/utau #Um*0.25/utau
alfaplus = params.L[1]/500.
betaplus = params.L[2]/200.
sigma = 0.00055 # 0.00055
epsilon = Um/200. #Um/200.
U[:] = 0
U[1] = Um*(Y-0.5*Y**2)
dev = 1+0.005*random.randn(Y.shape[0], Y.shape[1], Y.shape[2])
dd = utau*duplus/2.0*Xplus/40.*exp(-sigma*Xplus**2+0.5)*cos(betaplus*Zplus)*dev
U[1] += dd
U[2] += epsilon*sin(alfaplus*Yplus)*Xplus*exp(-sigma*Xplus**2)*dev
U_hat = solver.set_velocity(**context)
U = solver.get_velocity(**context)
U_hat = solver.set_velocity(**context)
if "KMM" in params.solver:
context.g[:] = 1j*context.K[1]*U_hat[2] - 1j*context.K[2]*U_hat[1]
# Set the flux
#flux[0] = context.FST.dx(U[1], context.ST.quad)
#solver.comm.Bcast(flux)
if solver.rank == 0:
print("Flux {}".format(flux[0]))
if not params.solver in ("KMM", "KMMRK3"):
P_hat = solver.compute_pressure(**context)
P = context.FST.ifst(P_hat, context.P, context.SN)
context.U_hat0[:] = context.U_hat[:]
context.H_hat1[:] = solver.get_convection(**context)
0
Example 78
Project: filterpy Source File: test_imm.py
def test_imm():
""" this test is drawn from Crassidis [1], example 4.6.
** References**
[1] Crassidis. "Optimal Estimation of Dynamic Systems", CRC Press,
Second edition.
"""
r = 1.
dt = 1.
phi_sim = np.array(
[[1, dt, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, dt],
[0, 0, 0, 1]])
gam = np.array([[dt**2/2, 0],
[dt, 0],
[0, dt**2/2],
[0, dt]])
x = np.array([[2000, 0, 10000, -15.]]).T
simxs = []
N = 600
for i in range(N):
x = np.dot(phi_sim, x)
if i >= 400:
x += np.dot(gam, np.array([[.075, .075]]).T)
simxs.append(x)
simxs = np.array(simxs)
#x = np.genfromtxt('c:/users/rlabbe/dropbox/Crassidis/mycode/x.csv', delimiter=',')
zs = np.zeros((N, 2))
for i in range(len(zs)):
zs[i, 0] = simxs[i, 0] + randn()*r
zs[i, 1] = simxs[i, 2] + randn()*r
try:
#data to test against crassidis' IMM matlab code
zs_tmp = np.genfromtxt('c:/users/rlabbe/dropbox/Crassidis/mycode/xx.csv', delimiter=',')[:-1]
zs = zs_tmp
except:
pass
ca = KalmanFilter(6, 2)
cano = KalmanFilter(6, 2)
dt2 = (dt**2)/2
ca.F = np.array(
[[1, dt, dt2, 0, 0, 0],
[0, 1, dt, 0, 0, 0],
[0, 0, 1, 0, 0, 0],
[0, 0, 0, 1, dt, dt2],
[0, 0, 0, 0, 1, dt],
[0, 0, 0, 0, 0, 1]])
cano.F = ca.F.copy()
ca.x = np.array([[2000., 0, 0, 10000, -15, 0]]).T
cano.x = ca.x.copy()
ca.P *= 1.e-12
cano.P *= 1.e-12
ca.R *= r**2
cano.R *= r**2
cano.Q *= 0
q = np.array([[.05, .125, 1/6],
[.125, 1/3, .5],
[1/6, .5, 1]])*1.e-3
ca.Q[0:3, 0:3] = q
ca.Q[3:6, 3:6] = q
ca.H = np.array([[1, 0, 0, 0, 0, 0],
[0, 0, 0, 1, 0, 0]])
cano.H = ca.H.copy()
filters = [ca, cano]
trans = np.array([[0.97, 0.03],
[0.03, 0.97]])
bank = IMMEstimator(filters, (0.5, 0.5), trans)
xs, probs = [], []
cvxs, caxs = [], []
for i, z in enumerate(zs[0:10]):
z = np.array([z]).T
bank.update(z)
#print(ca.likelihood, cano.likelihood)
#print(ca.x.T)
xs.append(bank.x.copy())
cvxs.append(ca.x.copy())
caxs.append(cano.x.copy())
#print(i, ca.likelihood, cano.likelihood, bank.w)
#print('p', bank.p)
probs.append(bank.mu.copy())
if DO_PLOT:
xs = np.array(xs)
cvxs = np.array(cvxs)
caxs = np.array(caxs)
probs = np.array(probs)
plt.subplot(121)
plt.plot(xs[:, 0], xs[:, 3], 'k')
#plt.plot(cvxs[:, 0], caxs[:, 3])
#plt.plot(simxs[:, 0], simxs[:, 2], 'g')
plt.scatter(zs[:, 0], zs[:, 1], marker='+', alpha=0.2)
plt.subplot(122)
plt.plot(probs[:, 0])
plt.plot(probs[:, 1])
plt.ylim(-1.5, 1.5)
plt.title('probability ratio p(cv)/p(ca)')
'''plt.figure()
0
Example 79
Project: pymanopt Source File: test_theano.py
def setUp(self):
x = T.vector()
y = T.matrix()
z = T.tensor3()
f = T.exp(T.sum(x**2)) + T.exp(T.sum(y**2)) + T.exp(T.sum(z**2))
self.cost = f
self.arg = [x, y, z]
n1 = self.n1 = 3
n2 = self.n2 = 4
n3 = self.n3 = 5
n4 = self.n4 = 6
n5 = self.n5 = 7
n6 = self.n6 = 8
self.y = y = (rnd.randn(n1), rnd.randn(n2, n3), rnd.randn(n4, n5, n6))
self.a = a = (rnd.randn(n1), rnd.randn(n2, n3), rnd.randn(n4, n5, n6))
self.correct_cost = (np.exp(np.sum(y[0]**2)) +
np.exp(np.sum(y[1]**2)) +
np.exp(np.sum(y[2]**2)))
# CALCULATE CORRECT GRAD
g1 = 2 * y[0] * np.exp(np.sum(y[0] ** 2))
g2 = 2 * y[1] * np.exp(np.sum(y[1] ** 2))
g3 = 2 * y[2] * np.exp(np.sum(y[2] ** 2))
self.correct_grad = (g1, g2, g3)
# CALCULATE CORRECT HESS
# 1. VECTOR
Ymat = np.matrix(y[0])
Amat = np.matrix(a[0])
diag = np.eye(n1)
H = np.exp(np.sum(y[0] ** 2)) * (4 * Ymat.T.dot(Ymat) + 2 * diag)
# Then 'left multiply' H by A
h1 = np.array(Amat.dot(H)).flatten()
# 2. MATRIX
# First form hessian tensor H (4th order)
Y1 = y[1].reshape(n2, n3, 1, 1)
Y2 = y[1].reshape(1, 1, n2, n3)
# Create an m x n x m x n array with diag[i,j,k,l] == 1 iff
# (i == k and j == l), this is a 'diagonal' tensor.
diag = np.eye(n2 * n3).reshape(n2, n3, n2, n3)
H = np.exp(np.sum(y[1] ** 2)) * (4 * Y1 * Y2 + 2 * diag)
# Then 'right multiply' H by A
Atensor = a[1].reshape(1, 1, n2, n3)
h2 = np.sum(H * Atensor, axis=(2, 3))
# 3. Tensor3
# First form hessian tensor H (6th order)
Y1 = y[2].reshape(n4, n5, n6, 1, 1, 1)
Y2 = y[2].reshape(1, 1, 1, n4, n5, n6)
# Create an n1 x n2 x n3 x n1 x n2 x n3 diagonal tensor
diag = np.eye(n4 * n5 * n6).reshape(n4, n5, n6, n4, n5, n6)
H = np.exp(np.sum(y[2] ** 2)) * (4 * Y1 * Y2 + 2 * diag)
# Then 'right multiply' H by A
Atensor = a[2].reshape(1, 1, 1, n4, n5, n6)
h3 = np.sum(H * Atensor, axis=(3, 4, 5))
self.correct_hess = (h1, h2, h3)
self.backend = TheanoBackend()
0
Example 80
Project: GPy Source File: dimensionality_reduction.py
def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=25, Q=4, sigma=.2):
import GPy
from pods.datasets import swiss_roll_generated
from GPy.models import BayesianGPLVM
data = swiss_roll_generated(num_samples=N, sigma=sigma)
Y = data['Y']
Y -= Y.mean()
Y /= Y.std()
t = data['t']
c = data['colors']
try:
from sklearn.manifold.isomap import Isomap
iso = Isomap().fit(Y)
X = iso.embedding_
if Q > 2:
X = _np.hstack((X, _np.random.randn(N, Q - 2)))
except ImportError:
X = _np.random.randn(N, Q)
if plot:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D # @UnusedImport
fig = plt.figure("Swiss Roll Data")
ax = fig.add_subplot(121, projection='3d')
ax.scatter(*Y.T, c=c)
ax.set_title("Swiss Roll")
ax = fig.add_subplot(122)
ax.scatter(*X.T[:2], c=c)
ax.set_title("BGPLVM init")
var = .5
S = (var * _np.ones_like(X) + _np.clip(_np.random.randn(N, Q) * var ** 2,
- (1 - var),
(1 - var))) + .001
Z = _np.random.permutation(X)[:num_inducing]
kernel = GPy.kern.RBF(Q, ARD=True) + GPy.kern.Bias(Q, _np.exp(-2)) + GPy.kern.White(Q, _np.exp(-2))
m = BayesianGPLVM(Y, Q, X=X, X_variance=S, num_inducing=num_inducing, Z=Z, kernel=kernel)
m.data_colors = c
m.data_t = t
if optimize:
m.optimize('bfgs', messages=verbose, max_iters=2e3)
if plot:
fig = plt.figure('fitted')
ax = fig.add_subplot(111)
s = m.input_sensitivity().argsort()[::-1][:2]
ax.scatter(*m.X.mean.T[s], c=c)
return m
0
Example 81
Project: GPy Source File: kernel_tests.py
def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verbose=False, fixed_X_dims=None):
"""
This function runs on kernels to check the correctness of their
implementation. It checks that the covariance function is positive definite
for a randomly generated data set.
:param kern: the kernel to be tested.
:type kern: GPy.kern.Kernpart
:param X: X input values to test the covariance function.
:type X: ndarray
:param X2: X2 input values to test the covariance function.
:type X2: ndarray
"""
pass_checks = True
if X is None:
X = np.random.randn(10, kern.input_dim)
if output_ind is not None:
X[:, output_ind] = np.random.randint(kern.output_dim, X.shape[0])
if X2 is None:
X2 = np.random.randn(20, kern.input_dim)
if output_ind is not None:
X2[:, output_ind] = np.random.randint(kern.output_dim, X2.shape[0])
if verbose:
print("Checking covariance function is positive definite.")
result = Kern_check_model(kern, X=X).is_positive_semi_definite()
if result and verbose:
print("Check passed.")
if not result:
print(("Positive definite check failed for " + kern.name + " covariance function."))
pass_checks = False
assert(result)
return False
if verbose:
print("Checking gradients of K(X, X) wrt theta.")
result = Kern_check_dK_dtheta(kern, X=X, X2=None).checkgrad(verbose=verbose)
if result and verbose:
print("Check passed.")
if not result:
print(("Gradient of K(X, X) wrt theta failed for " + kern.name + " covariance function. Gradient values as follows:"))
Kern_check_dK_dtheta(kern, X=X, X2=None).checkgrad(verbose=True)
pass_checks = False
assert(result)
return False
if verbose:
print("Checking gradients of K(X, X2) wrt theta.")
try:
result = Kern_check_dK_dtheta(kern, X=X, X2=X2).checkgrad(verbose=verbose)
except NotImplementedError:
result=True
if verbose:
print(("update_gradients_full, with differing X and X2, not implemented for " + kern.name))
if result and verbose:
print("Check passed.")
if not result:
print(("Gradient of K(X, X) wrt theta failed for " + kern.name + " covariance function. Gradient values as follows:"))
Kern_check_dK_dtheta(kern, X=X, X2=X2).checkgrad(verbose=True)
pass_checks = False
assert(result)
return False
if verbose:
print("Checking gradients of Kdiag(X) wrt theta.")
try:
result = Kern_check_dKdiag_dtheta(kern, X=X).checkgrad(verbose=verbose)
except NotImplementedError:
result=True
if verbose:
print(("update_gradients_diag not implemented for " + kern.name))
if result and verbose:
print("Check passed.")
if not result:
print(("Gradient of Kdiag(X) wrt theta failed for " + kern.name + " covariance function. Gradient values as follows:"))
Kern_check_dKdiag_dtheta(kern, X=X).checkgrad(verbose=True)
pass_checks = False
assert(result)
return False
if verbose:
print("Checking gradients of K(X, X) wrt X.")
try:
testmodel = Kern_check_dK_dX(kern, X=X, X2=None)
if fixed_X_dims is not None:
testmodel.X[:,fixed_X_dims].fix()
result = testmodel.checkgrad(verbose=verbose)
except NotImplementedError:
result=True
if verbose:
print(("gradients_X not implemented for " + kern.name))
if result and verbose:
print("Check passed.")
if not result:
print(("Gradient of K(X, X) wrt X failed for " + kern.name + " covariance function. Gradient values as follows:"))
testmodel.checkgrad(verbose=True)
assert(result)
pass_checks = False
return False
if verbose:
print("Checking gradients of K(X, X2) wrt X.")
try:
testmodel = Kern_check_dK_dX(kern, X=X, X2=X2)
if fixed_X_dims is not None:
testmodel.X[:,fixed_X_dims].fix()
result = testmodel.checkgrad(verbose=verbose)
except NotImplementedError:
result=True
if verbose:
print(("gradients_X not implemented for " + kern.name))
if result and verbose:
print("Check passed.")
if not result:
print(("Gradient of K(X, X2) wrt X failed for " + kern.name + " covariance function. Gradient values as follows:"))
testmodel.checkgrad(verbose=True)
assert(result)
pass_checks = False
return False
if verbose:
print("Checking gradients of Kdiag(X) wrt X.")
try:
testmodel = Kern_check_dKdiag_dX(kern, X=X)
if fixed_X_dims is not None:
testmodel.X[:,fixed_X_dims].fix()
result = testmodel.checkgrad(verbose=verbose)
except NotImplementedError:
result=True
if verbose:
print(("gradients_X not implemented for " + kern.name))
if result and verbose:
print("Check passed.")
if not result:
print(("Gradient of Kdiag(X) wrt X failed for " + kern.name + " covariance function. Gradient values as follows:"))
Kern_check_dKdiag_dX(kern, X=X).checkgrad(verbose=True)
pass_checks = False
assert(result)
return False
if verbose:
print("Checking gradients of dK(X, X2) wrt X2 with full cov in dimensions")
try:
testmodel = Kern_check_d2K_dXdX(kern, X=X, X2=X2)
if fixed_X_dims is not None:
testmodel.X[:,fixed_X_dims].fix()
result = testmodel.checkgrad(verbose=verbose)
except NotImplementedError:
result=True
if verbose:
print(("gradients_X not implemented for " + kern.name))
if result and verbose:
print("Check passed.")
if not result:
print(("Gradient of dK(X, X2) wrt X failed for " + kern.name + " covariance function. Gradient values as follows:"))
testmodel.checkgrad(verbose=True)
assert(result)
pass_checks = False
return False
if verbose:
print("Checking gradients of dK(X, X) wrt X with full cov in dimensions")
try:
testmodel = Kern_check_d2K_dXdX(kern, X=X, X2=None)
if fixed_X_dims is not None:
testmodel.X[:,fixed_X_dims].fix()
result = testmodel.checkgrad(verbose=verbose)
except NotImplementedError:
result=True
if verbose:
print(("gradients_X not implemented for " + kern.name))
if result and verbose:
print("Check passed.")
if not result:
print(("Gradient of dK(X, X) wrt X with full cov in dimensions failed for " + kern.name + " covariance function. Gradient values as follows:"))
testmodel.checkgrad(verbose=True)
assert(result)
pass_checks = False
return False
if verbose:
print("Checking gradients of dKdiag(X, X) wrt X with cov in dimensions")
try:
testmodel = Kern_check_d2Kdiag_dXdX(kern, X=X)
if fixed_X_dims is not None:
testmodel.X[:,fixed_X_dims].fix()
result = testmodel.checkgrad(verbose=verbose)
except NotImplementedError:
result=True
if verbose:
print(("gradients_X not implemented for " + kern.name))
if result and verbose:
print("Check passed.")
if not result:
print(("Gradient of dKdiag(X, X) wrt X with cov in dimensions failed for " + kern.name + " covariance function. Gradient values as follows:"))
testmodel.checkgrad(verbose=True)
assert(result)
pass_checks = False
return False
return pass_checks
0
Example 82
Project: ComputationalNeurodynamics Source File: Sync2Run.py
def Sync2Run():
"""
Docstring
"""
N1 = 800
N2 = 200
net = Sync2Connect(N1, N2)
T = 1000 # simulation time per episode
Ib = 5 # base current
# Initialise layers
for lr in xrange(net.Nlayers):
net.layer[lr].v = -65 * np.ones(net.layer[lr].N)
net.layer[lr].u = net.layer[lr].b * net.layer[lr].v
net.layer[lr].firings = np.array([])
# SIMULATE
for t in xrange(T):
# Deliver base current
net.layer[0].I = Ib*rn.randn(N1)
net.layer[1].I = Ib*rn.randn(N2)
# Delay the onset of activity for the second population
if t > 12:
net.layer[2].I = Ib*rn.randn(N1)
net.layer[3].I = Ib*rn.randn(N2)
else:
net.layer[2].I = np.zeros(N1)
net.layer[3].I = np.zeros(N2)
### Deliver a Poisson spike stream
# lambda = 0.01
# net.layer[1].I = 15*(poissrnd(lambda*15,N1) > 0)
# net.layer[2].I = 15*(poissrnd(lambda*15,N2) > 0)
# Update all the neurons
net.Update(t)
firings0 = net.layer[0].firings
firings2 = net.layer[2].firings
# Moving averages of firing rates in Hz for excitatory population
ws = 10 # window size
ds = 1 # slide window by ds
MF0 = np.zeros(np.ceil(T*1.0/ds))
MF2 = np.zeros(np.ceil(T*1.0/ds))
for j in xrange(1, T, ds):
MF0[j/ds] = sum(1*(firings0[:, 0] >= (j-ws)) * (firings0[:, 0] < j)) * 1000.0/(ws*N1)
MF2[j/ds] = sum(1*(firings2[:, 0] >= (j-ws)) * (firings2[:, 0] < j)) * 1000.0/(ws*N1)
# Raster plots of firings
plt.subplot(211)
if firings0.size is not 0:
plt.scatter(firings0[:, 0], firings0[:, 1], marker='.')
plt.xlim([0, T])
plt.ylabel('Neuron number')
plt.ylim([0, N1+1])
plt.subplot(212)
if firings2.size is not 0:
plt.scatter(firings2[:, 0], firings2[:, 1], marker='.')
plt.xlabel('Time (ms)')
plt.xlim([0, T])
plt.ylabel('Neuron number')
plt.ylim([0, N1+1])
# Plot mean firing rates
plt.figure()
plt.plot(range(0, T, ds), np.vstack([MF0, MF2]).T)
plt.xlabel('Time (ms)')
plt.ylabel('Mean firing rate')
plt.show()
return (MF0, MF2)
0
Example 83
Project: xarray Source File: test_indexing.py
def test_orthogonal_indexer(self):
x = np.random.randn(10, 11, 12, 13, 14)
y = np.arange(5)
I = ReturnItem()
# orthogonal and numpy indexing should be equivalent, because we only
# use at most one array and it never in between two slice objects
# (i.e., we try to avoid numpy's mind-boggling "partial indexing"
# http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html)
for i in [I[:], I[0], I[0, 0], I[:5], I[5:], I[2:5], I[3:-3], I[::-1],
I[::-2], I[5::-2], I[:3:-2], I[2:5:-1], I[7:3:-2], I[:3, :4],
I[:3, 0, :4], I[:3, 0, :4, 0], I[y], I[:, y], I[0, y],
I[:2, :3, y], I[0, y, :, :4, 0]]:
j = indexing.orthogonal_indexer(i, x.shape)
self.assertArrayEqual(x[i], x[j])
self.assertArrayEqual(self.set_to_zero(x, i),
self.set_to_zero(x, j))
# for more complicated cases, check orthogonal indexing is still
# equivalent to slicing
z = np.arange(2, 8, 2)
for i, j, shape in [
(I[y, y], I[:5, :5], (5, 5, 12, 13, 14)),
(I[y, z], I[:5, 2:8:2], (5, 3, 12, 13, 14)),
(I[0, y, y], I[0, :5, :5], (5, 5, 13, 14)),
(I[y, 0, z], I[:5, 0, 2:8:2], (5, 3, 13, 14)),
(I[y, :, z], I[:5, :, 2:8:2], (5, 11, 3, 13, 14)),
(I[0, :, z], I[0, :, 2:8:2], (11, 3, 13, 14)),
(I[0, :2, y, y, 0], I[0, :2, :5, :5, 0], (2, 5, 5)),
(I[0, :, y, :, 0], I[0, :, :5, :, 0], (11, 5, 13)),
(I[:, :, y, :, 0], I[:, :, :5, :, 0], (10, 11, 5, 13)),
(I[:, :, y, z, :], I[:, :, :5, 2:8:2], (10, 11, 5, 3, 14))]:
k = indexing.orthogonal_indexer(i, x.shape)
self.assertEqual(shape, x[k].shape)
self.assertArrayEqual(x[j], x[k])
self.assertArrayEqual(self.set_to_zero(x, j),
self.set_to_zero(x, k))
# standard numpy (non-orthogonal) indexing doesn't work anymore
with self.assertRaisesRegexp(ValueError, 'only supports 1d'):
indexing.orthogonal_indexer(x > 0, x.shape)
with self.assertRaisesRegexp(ValueError, 'invalid subkey'):
print(indexing.orthogonal_indexer((1.5 * y, 1.5 * y), x.shape))
0
Example 84
Project: python-control Source File: statesp.py
def _rss_generate(states, inputs, outputs, type):
"""Generate a random state space.
This does the actual random state space generation expected from rss and
drss. type is 'c' for continuous systems and 'd' for discrete systems.
"""
# Probability of repeating a previous root.
pRepeat = 0.05
# Probability of choosing a real root. Note that when choosing a complex
# root, the conjugate gets chosen as well. So the expected proportion of
# real roots is pReal / (pReal + 2 * (1 - pReal)).
pReal = 0.6
# Probability that an element in B or C will not be masked out.
pBCmask = 0.8
# Probability that an element in D will not be masked out.
pDmask = 0.3
# Probability that D = 0.
pDzero = 0.5
# Check for valid input arguments.
if states < 1 or states % 1:
raise ValueError("states must be a positive integer. states = %g." %
states)
if inputs < 1 or inputs % 1:
raise ValueError("inputs must be a positive integer. inputs = %g." %
inputs)
if outputs < 1 or outputs % 1:
raise ValueError("outputs must be a positive integer. outputs = %g." %
outputs)
# Make some poles for A. Preallocate a complex array.
poles = zeros(states) + zeros(states) * 0.j
i = 0
while i < states:
if rand() < pRepeat and i != 0 and i != states - 1:
# Small chance of copying poles, if we're not at the first or last
# element.
if poles[i-1].imag == 0:
# Copy previous real pole.
poles[i] = poles[i-1]
i += 1
else:
# Copy previous complex conjugate pair of poles.
poles[i:i+2] = poles[i-2:i]
i += 2
elif rand() < pReal or i == states - 1:
# No-oscillation pole.
if type == 'c':
poles[i] = -exp(randn()) + 0.j
elif type == 'd':
poles[i] = 2. * rand() - 1.
i += 1
else:
# Complex conjugate pair of oscillating poles.
if type == 'c':
poles[i] = complex(-exp(randn()), 3. * exp(randn()))
elif type == 'd':
mag = rand()
phase = 2. * pi * rand()
poles[i] = complex(mag * cos(phase),
mag * sin(phase))
poles[i+1] = complex(poles[i].real, -poles[i].imag)
i += 2
# Now put the poles in A as real blocks on the diagonal.
A = zeros((states, states))
i = 0
while i < states:
if poles[i].imag == 0:
A[i, i] = poles[i].real
i += 1
else:
A[i, i] = A[i+1, i+1] = poles[i].real
A[i, i+1] = poles[i].imag
A[i+1, i] = -poles[i].imag
i += 2
# Finally, apply a transformation so that A is not block-diagonal.
while True:
T = randn(states, states)
try:
A = dot(solve(T, A), T) # A = T \ A * T
break
except LinAlgError:
# In the unlikely event that T is rank-deficient, iterate again.
pass
# Make the remaining matrices.
B = randn(states, inputs)
C = randn(outputs, states)
D = randn(outputs, inputs)
# Make masks to zero out some of the elements.
while True:
Bmask = rand(states, inputs) < pBCmask
if any(Bmask): # Retry if we get all zeros.
break
while True:
Cmask = rand(outputs, states) < pBCmask
if any(Cmask): # Retry if we get all zeros.
break
if rand() < pDzero:
Dmask = zeros((outputs, inputs))
else:
Dmask = rand(outputs, inputs) < pDmask
# Apply masks.
B = B * Bmask
C = C * Cmask
D = D * Dmask
return StateSpace(A, B, C, D)
0
Example 85
def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tolerance=0.85, eps=1e-10, ax=None):
"""
Basic derivative check
Compares error decay of 0th and 1st order Taylor approximation at point
x0 for a randomized search direction.
:param callable fctn: function handle
:param numpy.array x0: point at which to check derivative
:param int num: number of times to reduce step length, h
:param bool plotIt: if you would like to plot
:param numpy.array dx: step direction
:param int expectedOrder: The order that you expect the derivative to yield.
:param float tolerance: The tolerance on the expected order.
:param float eps: What is zero?
:rtype: bool
:return: did you pass the test?!
.. plot::
:include-source:
from SimPEG import Tests, Utils, np
def simplePass(x):
return np.sin(x), Utils.sdiag(np.cos(x))
Tests.checkDerivative(simplePass, np.random.randn(5))
"""
print("{0!s} checkDerivative {1!s}".format('='*20, '='*20))
print("iter h |ft-f0| |ft-f0-h*J0*dx| Order\n{0!s}".format(('-'*57)))
f0, J0 = fctn(x0)
x0 = mkvc(x0)
if dx is None:
dx = np.random.randn(len(x0))
h = np.logspace(-1, -num, num)
E0 = np.ones(h.shape)
E1 = np.ones(h.shape)
def l2norm(x):
# because np.norm breaks if they are scalars?
return np.sqrt(np.real(np.vdot(x, x)))
for i in range(num):
# Evaluate at test point
ft, Jt = fctn( x0 + h[i]*dx )
# 0th order Taylor
E0[i] = l2norm( ft - f0 )
# 1st order Taylor
if inspect.isfunction(J0):
E1[i] = l2norm( ft - f0 - h[i]*J0(dx) )
else:
# We assume it is a numpy.ndarray
E1[i] = l2norm( ft - f0 - h[i]*J0.dot(dx) )
order0 = np.log10(E0[:-1]/E0[1:])
order1 = np.log10(E1[:-1]/E1[1:])
print(" {0:d} {1:1.2e} {2:1.3e} {3:1.3e} {4:1.3f}".format(i, h[i], E0[i], E1[i], np.nan if i == 0 else order1[i-1]))
# Ensure we are about precision
order0 = order0[E0[1:] > eps]
order1 = order1[E1[1:] > eps]
belowTol = (order1.size == 0 and order0.size >= 0)
# Make sure we get the correct order
correctOrder = order1.size > 0 and np.mean(order1) > tolerance * expectedOrder
passTest = belowTol or correctOrder
if passTest:
print("{0!s} PASS! {1!s}".format('='*25, '='*25))
print(happiness[np.random.randint(len(happiness))]+'\n')
else:
print("{0!s}\n{1!s} FAIL! {2!s}\n{3!s}".format('*'*57, '<'*25, '>'*25, '*'*57))
print(sadness[np.random.randint(len(sadness))]+'\n')
if plotIt:
import matplotlib.pyplot as plt
ax = ax or plt.subplot(111)
ax.loglog(h, E0, 'b')
ax.loglog(h, E1, 'g--')
ax.set_title('Check Derivative - {0!s}'.format(('PASSED :)' if passTest else 'FAILED :(')))
ax.set_xlabel('h')
ax.set_ylabel('Error')
leg = ax.legend(['$\mathcal{O}(h)$', '$\mathcal{O}(h^2)$'], loc='best',
title="$f(x + h\Delta x) - f(x) - h g(x) \Delta x - \mathcal{O}(h^2) = 0$",
frameon=False)
plt.setp(leg.get_title(),fontsize=15)
plt.show()
return passTest
0
Example 86
Project: theano_pyglm Source File: generate_synth_data.py
def run_gen_synth_data():
""" Run a test with synthetic data and MCMC inference
"""
options, args = parse_cmd_line_args()
# Create the model
dt = 0.001
model = make_model(options.model, N=options.N, dt=dt)
# Set the sparsity level to minimize the risk of unstable networks
stabilize_sparsity(model)
print "Creating master population object"
popn = Population(model)
# Sample random parameters from the model
x_true = popn.sample()
# Check stability of matrix
assert check_stability(model, x_true, options.N), "ERROR: Sampled network is unstable!"
# Save the model so it can be loaded alongside the data
fname_model = os.path.join(options.resultsDir, 'model.pkl')
print "Saving data to %s" % fname_model
with open(fname_model,'w') as f:
cPickle.dump(model, f, protocol=-1)
print "Generating synthetic data with %d neurons and %.2f seconds." % \
(options.N, options.T_stop)
# Set simulation parametrs
dt_stim = 0.1
D_stim = (5,5)
# D_stim = model['bkgd']['D_stim'] if 'D_stim' in model['bkgd'] else 0
if isinstance(D_stim, int):
D_stim = [D_stim]
stim = np.random.randn(options.T_stop/dt_stim, *D_stim)
data = gen_synth_data(options.N, options.T_stop, popn, x_true, dt, dt_stim, D_stim, stim)
# Set the data so that the population state can be evaluated
popn.add_data(data)
# DEBUG Evaluate the firing rate and the simulated firing rate
state = popn.eval_state(x_true)
for n in np.arange(options.N):
lam_true = state['glms'][n]['lam']
lam_sim = popn.glm.nlin_model.f_nlin(data['X'][:,n])
assert np.allclose(lam_true, lam_sim)
# Pickle the data so we can open it more easily
fname_pkl = os.path.join(options.resultsDir, 'data.pkl')
print "Saving data to %s" % fname_pkl
with open(fname_pkl,'w') as f:
cPickle.dump(data, f, protocol=-1)
# Plot firing rates, stimulus responses, etc
do_plot_imp_resonses = int(options.N) <= 16
plot_results(popn, data['vars'],
resdir=options.resultsDir,
do_plot_stim_resp=True,
do_plot_imp_responses=do_plot_imp_resonses)
0
Example 87
def simulate(self, ts_length=100):
"""
Simulate a time series of length ts_length, first drawing
x_0 ~ N(mu_0, Sigma_0)
Parameters
----------
ts_length : scalar(int), optional(default=100)
The length of the simulation
Returns
-------
x : array_like(float)
An n x ts_length array, where the t-th column is x_t
y : array_like(float)
A k x ts_length array, where the t-th column is y_t
"""
x0 = multivariate_normal(self.mu_0.flatten(), self.Sigma_0)
w = np.random.randn(self.m, ts_length-1)
v = self.C.dot(w) # Multiply each w_t by C to get v_t = C w_t
# == simulate time series == #
x = simulate_linear_model(self.A, x0, v, ts_length)
if self.H is not None:
v = np.random.randn(self.l, ts_length)
y = self.G.dot(x) + self.H.dot(v)
else:
y = self.G.dot(x)
return x, y
0
Example 88
Project: bayespy Source File: test_dot.py
def test_message_to_child(self):
"""
Test the message from SumMultiply to its children.
"""
def compare_moments(u0, u1, *args):
Y = SumMultiply(*args)
u_Y = Y.get_moments()
self.assertAllClose(u_Y[0], u0)
self.assertAllClose(u_Y[1], u1)
# Test constant parent
y = np.random.randn(2,3,4)
compare_moments(y,
linalg.outer(y, y, ndim=2),
'ij->ij',
y)
# Do nothing for 2-D array
Y = GaussianARD(np.random.randn(5,2,3),
np.random.rand(5,2,3),
plates=(5,),
shape=(2,3))
y = Y.get_moments()
compare_moments(y[0],
y[1],
'ij->ij',
Y)
compare_moments(y[0],
y[1],
Y,
[0,1],
[0,1])
# Sum over the rows of a matrix
Y = GaussianARD(np.random.randn(5,2,3),
np.random.rand(5,2,3),
plates=(5,),
shape=(2,3))
y = Y.get_moments()
mu = np.einsum('...ij->...j', y[0])
cov = np.einsum('...ijkl->...jl', y[1])
compare_moments(mu,
cov,
'ij->j',
Y)
compare_moments(mu,
cov,
Y,
[0,1],
[1])
# Inner product of three vectors
X1 = GaussianARD(np.random.randn(2),
np.random.rand(2),
plates=(),
shape=(2,))
x1 = X1.get_moments()
X2 = GaussianARD(np.random.randn(6,1,2),
np.random.rand(6,1,2),
plates=(6,1),
shape=(2,))
x2 = X2.get_moments()
X3 = GaussianARD(np.random.randn(7,6,5,2),
np.random.rand(7,6,5,2),
plates=(7,6,5),
shape=(2,))
x3 = X3.get_moments()
mu = np.einsum('...i,...i,...i->...', x1[0], x2[0], x3[0])
cov = np.einsum('...ij,...ij,...ij->...', x1[1], x2[1], x3[1])
compare_moments(mu,
cov,
'i,i,i',
X1,
X2,
X3)
compare_moments(mu,
cov,
'i,i,i->',
X1,
X2,
X3)
compare_moments(mu,
cov,
X1,
[9],
X2,
[9],
X3,
[9])
compare_moments(mu,
cov,
X1,
[9],
X2,
[9],
X3,
[9],
[])
# Outer product of two vectors
X1 = GaussianARD(np.random.randn(2),
np.random.rand(2),
plates=(5,),
shape=(2,))
x1 = X1.get_moments()
X2 = GaussianARD(np.random.randn(6,1,2),
np.random.rand(6,1,2),
plates=(6,1),
shape=(2,))
x2 = X2.get_moments()
mu = np.einsum('...i,...j->...ij', x1[0], x2[0])
cov = np.einsum('...ik,...jl->...ijkl', x1[1], x2[1])
compare_moments(mu,
cov,
'i,j->ij',
X1,
X2)
compare_moments(mu,
cov,
X1,
[9],
X2,
[7],
[9,7])
# Matrix product
Y1 = GaussianARD(np.random.randn(3,2),
np.random.rand(3,2),
plates=(),
shape=(3,2))
y1 = Y1.get_moments()
Y2 = GaussianARD(np.random.randn(5,2,3),
np.random.rand(5,2,3),
plates=(5,),
shape=(2,3))
y2 = Y2.get_moments()
mu = np.einsum('...ik,...kj->...ij', y1[0], y2[0])
cov = np.einsum('...ikjl,...kmln->...imjn', y1[1], y2[1])
compare_moments(mu,
cov,
'ik,kj->ij',
Y1,
Y2)
compare_moments(mu,
cov,
Y1,
['i','k'],
Y2,
['k','j'],
['i','j'])
# Trace of a matrix product
Y1 = GaussianARD(np.random.randn(3,2),
np.random.rand(3,2),
plates=(),
shape=(3,2))
y1 = Y1.get_moments()
Y2 = GaussianARD(np.random.randn(5,2,3),
np.random.rand(5,2,3),
plates=(5,),
shape=(2,3))
y2 = Y2.get_moments()
mu = np.einsum('...ij,...ji->...', y1[0], y2[0])
cov = np.einsum('...ikjl,...kilj->...', y1[1], y2[1])
compare_moments(mu,
cov,
'ij,ji',
Y1,
Y2)
compare_moments(mu,
cov,
'ij,ji->',
Y1,
Y2)
compare_moments(mu,
cov,
Y1,
['i','j'],
Y2,
['j','i'])
compare_moments(mu,
cov,
Y1,
['i','j'],
Y2,
['j','i'],
[])
# Vector-matrix-vector product
X1 = GaussianARD(np.random.randn(3),
np.random.rand(3),
plates=(),
shape=(3,))
x1 = X1.get_moments()
X2 = GaussianARD(np.random.randn(6,1,2),
np.random.rand(6,1,2),
plates=(6,1),
shape=(2,))
x2 = X2.get_moments()
Y = GaussianARD(np.random.randn(3,2),
np.random.rand(3,2),
plates=(),
shape=(3,2))
y = Y.get_moments()
mu = np.einsum('...i,...ij,...j->...', x1[0], y[0], x2[0])
cov = np.einsum('...ia,...ijab,...jb->...', x1[1], y[1], x2[1])
compare_moments(mu,
cov,
'i,ij,j',
X1,
Y,
X2)
compare_moments(mu,
cov,
X1,
[1],
Y,
[1,2],
X2,
[2])
# Complex sum-product of 0-D, 1-D, 2-D and 3-D arrays
V = GaussianARD(np.random.randn(7,6,5),
np.random.rand(7,6,5),
plates=(7,6,5),
shape=())
v = V.get_moments()
X = GaussianARD(np.random.randn(6,1,2),
np.random.rand(6,1,2),
plates=(6,1),
shape=(2,))
x = X.get_moments()
Y = GaussianARD(np.random.randn(3,4),
np.random.rand(3,4),
plates=(5,),
shape=(3,4))
y = Y.get_moments()
Z = GaussianARD(np.random.randn(4,2,3),
np.random.rand(4,2,3),
plates=(6,5),
shape=(4,2,3))
z = Z.get_moments()
mu = np.einsum('...,...i,...kj,...jik->...k', v[0], x[0], y[0], z[0])
cov = np.einsum('...,...ia,...kjcb,...jikbac->...kc', v[1], x[1], y[1], z[1])
compare_moments(mu,
cov,
',i,kj,jik->k',
V,
X,
Y,
Z)
compare_moments(mu,
cov,
V,
[],
X,
['i'],
Y,
['k','j'],
Z,
['j','i','k'],
['k'])
#
# Gaussian-gamma parents
#
# Outer product of vectors
X1 = GaussianARD(np.random.randn(2),
np.random.rand(2),
shape=(2,))
x1 = X1.get_moments()
X2 = GaussianGamma(
np.random.randn(6,1,2),
random.covariance(2),
np.random.rand(6,1),
np.random.rand(6,1),
plates=(6,1)
)
x2 = X2.get_moments()
Y = SumMultiply('i,j->ij', X1, X2)
u = Y._message_to_child()
y = np.einsum('...i,...j->...ij', x1[0], x2[0])
yy = np.einsum('...ik,...jl->...ijkl', x1[1], x2[1])
self.assertAllClose(u[0], y)
self.assertAllClose(u[1], yy)
self.assertAllClose(u[2], x2[2])
self.assertAllClose(u[3], x2[3])
pass
0
Example 89
Project: bayespy Source File: test_gaussian_markov_chain.py
def test_message_to_parents_with_inputs(self):
""" Check gradient passed to inputs parent node """
def check(Mu, Lambda, A, V, U):
X = GaussianMarkovChain(Mu, Lambda, A, V, inputs=U)
Y = Gaussian(X, random.covariance(D))
# Check moments
self.assert_moments(
X,
postprocess=lambda u: [
u[0],
u[1] + linalg.transpose(u[1], ndim=1),
u[2]
]
)
Y.observe(np.random.randn(N+1, D))
X.update()
# Check gradient messages to parents
self.assert_message_to_parent(X, Mu)
self.assert_message_to_parent(
X,
Lambda,
postprocess=lambda phi: [
phi[0] + linalg.transpose(phi[0], ndim=1),
phi[1]
]
)
self.assert_message_to_parent(
X,
A,
postprocess=lambda phi: [
phi[0],
phi[1] + linalg.transpose(phi[1], ndim=1),
]
)
self.assert_message_to_parent(X, V)
self.assert_message_to_parent(X, U)
N = 4
D = 2
K = 3
check(
Gaussian(
np.random.randn(D),
random.covariance(D)
),
Wishart(
D,
random.covariance(D)
),
Gaussian(
np.random.randn(D,D+K),
random.covariance(D+K)
),
Gamma(
D,
np.random.rand(D)
),
Gaussian(
np.random.randn(N,K),
random.covariance(K)
)
)
check(
Gaussian(
np.random.randn(D),
random.covariance(D)
),
Wishart(
D,
random.covariance(D)
),
GaussianGamma(
np.random.randn(D,D+K),
random.covariance(D+K),
D,
np.random.rand(D),
ndim=1
),
Gamma(
D,
np.random.rand(D)
),
Gaussian(
np.random.randn(N,K),
random.covariance(K)
)
)
pass
0
Example 90
Project: zipline Source File: test_factor.py
@parameterized.expand(gen_ranking_cases())
def test_masked_rankdata_2d(self,
seed_value,
method,
use_mask,
set_missing,
ascending):
eyemask = ~eye(5, dtype=bool)
nomask = ones((5, 5), dtype=bool)
seed(seed_value)
asfloat = (randn(5, 5) * seed_value)
asdatetime = (asfloat).copy().view('datetime64[ns]')
mask = eyemask if use_mask else nomask
if set_missing:
asfloat[:, 2] = nan
asdatetime[:, 2] = NaTns
float_result = masked_rankdata_2d(
data=asfloat,
mask=mask,
missing_value=nan,
method=method,
ascending=True,
)
datetime_result = masked_rankdata_2d(
data=asdatetime,
mask=mask,
missing_value=NaTns,
method=method,
ascending=True,
)
check_arrays(float_result, datetime_result)
0
Example 91
def test_message_to_child(self):
"""
Test the updating of GaussianMarkovChain.
Check that the moments and the lower bound contribution are computed
correctly.
"""
# TODO: Add plates and missing values!
# Dimensionalities
D = 3
N = 5
(Y, X, Mu, Lambda, A, V) = self.create_model(N, D)
# Inference with arbitrary observations
y = np.random.randn(N,D)
Y.observe(y)
X.update()
(x_vb, xnxn_vb, xpxn_vb) = X.get_moments()
# Get parameter moments
(mu0, mumu0) = Mu.get_moments()
(icov0, logdet0) = Lambda.get_moments()
(a, aa) = A.get_moments()
(icov_x, logdetx) = V.get_moments()
icov_x = np.diag(icov_x)
# Prior precision
Z = np.einsum('...kij,...kk->...ij', aa, icov_x)
U_diag = [icov0+Z] + (N-2)*[icov_x+Z] + [icov_x]
U_super = (N-1) * [-np.dot(a.T, icov_x)]
U = misc.block_banded(U_diag, U_super)
# Prior mean
mu_prior = np.zeros(D*N)
mu_prior[:D] = np.dot(icov0,mu0)
# Data
Cov = np.linalg.inv(U + np.identity(D*N))
mu = np.dot(Cov, mu_prior + y.flatten())
# Moments
xx = mu[:,np.newaxis]*mu[np.newaxis,:] + Cov
mu = np.reshape(mu, (N,D))
xx = np.reshape(xx, (N,D,N,D))
# Check results
self.assertAllClose(x_vb, mu,
msg="Incorrect mean")
for n in range(N):
self.assertAllClose(xnxn_vb[n,:,:], xx[n,:,n,:],
msg="Incorrect second moment")
for n in range(N-1):
self.assertAllClose(xpxn_vb[n,:,:], xx[n,:,n+1,:],
msg="Incorrect lagged second moment")
# Compute the entropy H(X)
ldet = linalg.logdet_cov(Cov)
H = random.gaussian_entropy(-ldet, N*D)
# Compute <log p(X|...)>
xx = np.reshape(xx, (N*D, N*D))
mu = np.reshape(mu, (N*D,))
ldet = -logdet0 - np.sum(np.ones((N-1,D))*logdetx)
P = random.gaussian_logpdf(np.einsum('...ij,...ij',
xx,
U),
np.einsum('...i,...i',
mu,
mu_prior),
np.einsum('...ij,...ij',
mumu0,
icov0),
-ldet,
N*D)
# The VB bound from the net
l = X.lower_bound_contribution()
self.assertAllClose(l, H+P)
# Compute the true bound <log p(X|...)> + H(X)
#
# Simple tests
#
def check(N, D, plates=None, mu=None, Lambda=None, A=None, V=None):
if mu is None:
mu = np.random.randn(D)
if Lambda is None:
Lambda = random.covariance(D)
if A is None:
A = np.random.randn(D,D)
if V is None:
V = np.random.rand(D)
X = GaussianMarkovChain(mu,
Lambda,
A,
V,
plates=plates,
n=N)
(u0, u1, u2) = X._message_to_child()
(mu, mumu) = Gaussian._ensure_moments(mu, GaussianMoments, ndim=1).get_moments()
(Lambda, _) = Wishart._ensure_moments(Lambda, WishartMoments, ndim=1).get_moments()
(a, aa) = Gaussian._ensure_moments(A, GaussianMoments, ndim=1).get_moments()
a = a * np.ones((N-1,D,D)) # explicit broadcasting for simplicity
aa = aa * np.ones((N-1,D,D,D)) # explicit broadcasting for simplicity
(v, _) = Gamma._ensure_moments(V, GammaMoments).get_moments()
v = v * np.ones((N-1,D))
plates_C = X.plates
plates_mu = X.plates
C = np.zeros(plates_C + (N,D,N,D))
plates_mu = np.shape(mu)[:-1]
m = np.zeros(plates_mu + (N,D))
m[...,0,:] = np.einsum('...ij,...j->...i', Lambda, mu)
C[...,0,:,0,:] = Lambda + np.einsum('...dij,...d->...ij',
aa[...,0,:,:,:],
v[...,0,:])
for n in range(1,N-1):
C[...,n,:,n,:] = (np.einsum('...dij,...d->...ij',
aa[...,n,:,:,:],
v[...,n,:])
+ v[...,n,:,None] * np.identity(D))
for n in range(N-1):
C[...,n,:,n+1,:] = -np.einsum('...di,...d->...id',
a[...,n,:,:],
v[...,n,:])
C[...,n+1,:,n,:] = -np.einsum('...di,...d->...di',
a[...,n,:,:],
v[...,n,:])
C[...,-1,:,-1,:] = v[...,-1,:,None]*np.identity(D)
C = np.reshape(C, plates_C+(N*D,N*D))
Cov = np.linalg.inv(C)
Cov = np.reshape(Cov, plates_C+(N,D,N,D))
m0 = np.einsum('...minj,...nj->...mi', Cov, m)
m1 = np.zeros(plates_C+(N,D,D))
m2 = np.zeros(plates_C+(N-1,D,D))
for n in range(N):
m1[...,n,:,:] = Cov[...,n,:,n,:] + np.einsum('...i,...j->...ij',
m0[...,n,:],
m0[...,n,:])
for n in range(N-1):
m2[...,n,:,:] = Cov[...,n,:,n+1,:] + np.einsum('...i,...j->...ij',
m0[...,n,:],
m0[...,n+1,:])
self.assertAllClose(m0, u0*np.ones(np.shape(m0)))
self.assertAllClose(m1, u1*np.ones(np.shape(m1)))
self.assertAllClose(m2, u2*np.ones(np.shape(m2)))
pass
check(4,1)
check(4,3)
#
# Test mu
#
# Simple
check(4,3,
mu=Gaussian(np.random.randn(3),
random.covariance(3)))
# Plates
check(4,3,
mu=Gaussian(np.random.randn(5,6,3),
random.covariance(3),
plates=(5,6)))
# Plates with moments broadcasted over plates
check(4,3,
mu=Gaussian(np.random.randn(3),
random.covariance(3),
plates=(5,)))
check(4,3,
mu=Gaussian(np.random.randn(1,3),
random.covariance(3),
plates=(5,)))
# Plates broadcasting
check(4,3,
plates=(5,),
mu=Gaussian(np.random.randn(3),
random.covariance(3),
plates=()))
check(4,3,
plates=(5,),
mu=Gaussian(np.random.randn(1,3),
random.covariance(3),
plates=(1,)))
#
# Test Lambda
#
# Simple
check(4,3,
Lambda=Wishart(10+np.random.rand(),
random.covariance(3)))
# Plates
check(4,3,
Lambda=Wishart(10+np.random.rand(),
random.covariance(3),
plates=(5,6)))
# Plates with moments broadcasted over plates
check(4,3,
Lambda=Wishart(10+np.random.rand(),
random.covariance(3),
plates=(5,)))
check(4,3,
Lambda=Wishart(10+np.random.rand(1),
random.covariance(3),
plates=(5,)))
# Plates broadcasting
check(4,3,
plates=(5,),
Lambda=Wishart(10+np.random.rand(),
random.covariance(3),
plates=()))
check(4,3,
plates=(5,),
Lambda=Wishart(10+np.random.rand(),
random.covariance(3),
plates=(1,)))
#
# Test A
#
# Simple
check(4,3,
A=GaussianARD(np.random.randn(3,3),
np.random.rand(3,3),
shape=(3,),
plates=(3,)))
# Plates on time axis
check(5,3,
A=GaussianARD(np.random.randn(4,3,3),
np.random.rand(4,3,3),
shape=(3,),
plates=(4,3)))
# Plates on time axis with broadcasted moments
check(5,3,
A=GaussianARD(np.random.randn(1,3,3),
np.random.rand(1,3,3),
shape=(3,),
plates=(4,3)))
check(5,3,
A=GaussianARD(np.random.randn(3,3),
np.random.rand(3,3),
shape=(3,),
plates=(4,3)))
# Plates
check(4,3,
A=GaussianARD(np.random.randn(5,6,1,3,3),
np.random.rand(5,6,1,3,3),
shape=(3,),
plates=(5,6,1,3)))
# Plates with moments broadcasted over plates
check(4,3,
A=GaussianARD(np.random.randn(3,3),
np.random.rand(3,3),
shape=(3,),
plates=(5,1,3)))
check(4,3,
A=GaussianARD(np.random.randn(1,1,3,3),
np.random.rand(1,1,3,3),
shape=(3,),
plates=(5,1,3)))
# Plates broadcasting
check(4,3,
plates=(5,),
A=GaussianARD(np.random.randn(3,3),
np.random.rand(3,3),
shape=(3,),
plates=(3,)))
check(4,3,
plates=(5,),
A=GaussianARD(np.random.randn(3,3),
np.random.rand(3,3),
shape=(3,),
plates=(1,1,3)))
#
# Test v
#
# Simple
check(4,3,
V=Gamma(np.random.rand(1,3),
np.random.rand(1,3),
plates=(1,3)))
check(4,3,
V=Gamma(np.random.rand(3),
np.random.rand(3),
plates=(3,)))
# Plates
check(4,3,
V=Gamma(np.random.rand(5,6,1,3),
np.random.rand(5,6,1,3),
plates=(5,6,1,3)))
# Plates with moments broadcasted over plates
check(4,3,
V=Gamma(np.random.rand(1,3),
np.random.rand(1,3),
plates=(5,1,3)))
check(4,3,
V=Gamma(np.random.rand(1,1,3),
np.random.rand(1,1,3),
plates=(5,1,3)))
# Plates broadcasting
check(4,3,
plates=(5,),
V=Gamma(np.random.rand(1,3),
np.random.rand(1,3),
plates=(1,3)))
check(4,3,
plates=(5,),
V=Gamma(np.random.rand(1,1,3),
np.random.rand(1,1,3),
plates=(1,1,3)))
#
# Check with input signals
#
mu = 2
Lambda = 3
A = 4
B = 5
v = 6
inputs = [[-2], [3]]
X = GaussianMarkovChain([mu], [[Lambda]], [[A,B]], [v], inputs=inputs)
V = (np.array([[v*A**2, -v*A, 0],
[-v*A, v*A**2, -v*A],
[0, -v*A, 0]]) +
np.array([[Lambda, 0, 0],
[0, v, 0],
[0, 0, v]]))
m = (np.array([Lambda*mu, 0, 0]) +
np.array([0, v*B*inputs[0][0], v*B*inputs[1][0]]) -
np.array([v*A*B*inputs[0][0], v*A*B*inputs[1][0], 0]))
Cov = np.linalg.inv(V)
mean = np.dot(Cov, m)
X.update()
u = X.get_moments()
self.assertAllClose(u[0], mean[:,None])
self.assertAllClose(u[1] - u[0][...,None,:]*u[0][...,:,None],
Cov[(0,1,2),(0,1,2),None,None])
self.assertAllClose(u[2] - u[0][...,:-1,:,None]*u[0][...,1:,None,:],
Cov[(0,1),(1,2),None,None])
pass
0
Example 92
Project: bayespy Source File: test_node.py
def test_compute_message(self):
"""
Test the general sum-multiply function for message computations
"""
self.assertAllClose(Node._compute_message(3,
plates_from=(),
plates_to=(),
ndim=0),
3)
# Sum over one array
self.assertAllClose(Node._compute_message([1, 2, 3],
plates_from=(3,),
plates_to=(),
ndim=0),
6)
# Sum plates
self.assertAllClose(Node._compute_message([1, 2, 3],
[4, 4, 4],
[5, 5, 5],
plates_from=(3,),
plates_to=(),
ndim=0),
20+40+60)
# Do not sum plates
self.assertAllClose(Node._compute_message([1, 2, 3],
[4, 4, 4],
[5, 5, 5],
plates_from=(3,),
plates_to=(3,),
ndim=0),
[20, 40, 60])
# Give ndim
self.assertAllClose(Node._compute_message([1, 2, 3],
[4, 4, 4],
[5, 5, 5],
plates_from=(),
plates_to=(),
ndim=1),
[20, 40, 60])
# Broadcast plates_from
self.assertAllClose(Node._compute_message(3,
4,
5,
plates_from=(3,),
plates_to=(),
ndim=0),
3 * (3*4*5))
# Broadcast plates_to
self.assertAllClose(Node._compute_message(3,
4,
5,
plates_from=(3,),
plates_to=(3,),
ndim=0),
3*4*5)
# Different ndims
self.assertAllClose(Node._compute_message([1, 2, 3],
[4],
5,
plates_from=(3,),
plates_to=(3,),
ndim=0),
[1*4*5, 2*4*5, 3*4*5])
# Broadcasting dims for some arrays
self.assertAllClose(Node._compute_message([1, 2, 3],
[4],
5,
plates_from=(),
plates_to=(),
ndim=1),
[1*4*5, 2*4*5, 3*4*5])
# Bugfix: Check that plate keys are mapped correctly
self.assertAllClose(
Node._compute_message(
[[1], [2], [3]],
plates_from=(3,2),
plates_to=(1,2),
ndim=0
),
[[6]]
)
# Bugfix: Check plate key mapping when plates_to is shorter than shape
# of the array
self.assertAllClose(
Node._compute_message(
[[1, 2, 3], [4, 5, 6]],
plates_from=(2,3),
plates_to=(3,),
ndim=0
),
[5, 7, 9]
)
# Complex example
x1 = np.random.randn(5,4,1,2,1)
x2 = np.random.randn( 1,2,1)
x3 = np.random.randn(5,1,1,1,1)
self.assertAllClose(Node._compute_message(x1, x2, x3,
plates_from=(6,5,4,3),
plates_to=(5,1,1),
ndim=2),
3*6*np.sum(x1*x2*x3, axis=(-4,-3), keepdims=True))
pass
0
Example 93
def tsne(X = Math.array([]), no_dims = 2, initial_dims = 50, perplexity = 30.0):
"""Runs t-SNE on the dataset in the NxD array X to reduce its dimensionality to no_dims dimensions.
The syntaxis of the function is Y = tsne.tsne(X, no_dims, perplexity), where X is an NxD NumPy array."""
# Check inputs
if X.dtype != "float64":
print "Error: array X should have type float64.";
return -1;
#if no_dims.__class__ != "<type 'int'>": # doesn't work yet!
# print "Error: number of dimensions should be an integer.";
# return -1;
# Initialize variables
X = pca(X, initial_dims).real;
(n, d) = X.shape;
max_iter = 1000;
initial_momentum = 0.5;
final_momentum = 0.8;
eta = 500;
min_gain = 0.01;
Y = Math.random.randn(n, no_dims);
dY = Math.zeros((n, no_dims));
iY = Math.zeros((n, no_dims));
gains = Math.ones((n, no_dims));
# Compute P-values
P = x2p(X, 1e-5, perplexity);
P = P + Math.transpose(P);
P = P / Math.sum(P);
P = P * 4; # early exaggeration
P = Math.maximum(P, 1e-12);
# Run iterations
for iter in range(max_iter):
# Compute pairwise affinities
sum_Y = Math.sum(Math.square(Y), 1);
num = 1 / (1 + Math.add(Math.add(-2 * Math.dot(Y, Y.T), sum_Y).T, sum_Y));
num[range(n), range(n)] = 0;
Q = num / Math.sum(num);
Q = Math.maximum(Q, 1e-12);
# Compute gradient
PQ = P - Q;
for i in range(n):
dY[i,:] = Math.sum(Math.tile(PQ[:,i] * num[:,i], (no_dims, 1)).T * (Y[i,:] - Y), 0);
# Perform the update
if iter < 20:
momentum = initial_momentum
else:
momentum = final_momentum
gains = (gains + 0.2) * ((dY > 0) != (iY > 0)) + (gains * 0.8) * ((dY > 0) == (iY > 0));
gains[gains < min_gain] = min_gain;
iY = momentum * iY - eta * (gains * dY);
Y = Y + iY;
Y = Y - Math.tile(Math.mean(Y, 0), (n, 1));
# Compute current value of cost function
if (iter + 1) % 10 == 0:
C = Math.sum(P * Math.log(P / Q));
print "Iteration ", (iter + 1), ": error is ", C
# Stop lying about P-values
if iter == 100:
P = P / 4;
# Return solution
return Y;
0
Example 94
def test_message_to_child(self):
"""
Test message to child of X[..] node operator.
"""
class DummyNode(Node):
_moments = Moments()
_parent_moments = (Moments(),)
def __init__(self, u, **kwargs):
self.u = u
super().__init__(**kwargs)
def _message_to_child(self):
return self.u
def _get_id_list(self):
return []
# Message not a reference to X.u but a copy of it
X = DummyNode([np.random.randn(3)],
plates=(3,),
dims=((),))
Y = X[2]
self.assertTrue(Y._message_to_child() is not X.u,
msg="Slice node operator sends a reference to the "
"node's moment list as a message instead of a copy "
"of the list.")
# Integer indices
X = DummyNode([np.random.randn(3,4)],
plates=(3,4),
dims=((),))
Y = X[2,1]
self.assertMessageToChild(Y, [X.u[0][2,1]])
# Too few integer indices
X = DummyNode([np.random.randn(3,4)],
plates=(3,4),
dims=((),))
Y = X[2]
self.assertMessageToChild(Y, [X.u[0][2]])
# Integer for broadcasted moment
X = DummyNode([np.random.randn(4)],
plates=(3,4),
dims=((),))
Y = X[2,1]
self.assertMessageToChild(Y, [X.u[0][1]])
X = DummyNode([np.random.randn(4,1)],
plates=(3,4),
dims=((),))
Y = X[2,1]
self.assertMessageToChild(Y, [X.u[0][2,0]])
# Ignore leading new axes
X = DummyNode([np.random.randn(3)],
plates=(3,),
dims=((),))
Y = X[None,None,2]
self.assertMessageToChild(Y, [X.u[0][2]])
# Ignore new axes before missing+broadcasted plate axes
X = DummyNode([np.random.randn(3)],
plates=(4,3,),
dims=((),))
Y = X[1,None,None,2]
self.assertMessageToChild(Y, [X.u[0][2]])
# New axes
X = DummyNode([np.random.randn(3,4)],
plates=(3,4),
dims=((),))
Y = X[2,None,None,1]
self.assertMessageToChild(Y, [X.u[0][2,None,None,1]])
# New axes for broadcasted axes
X = DummyNode([np.random.randn(4)],
plates=(3,4),
dims=((),))
Y = X[2,1,None,None]
self.assertMessageToChild(Y, [X.u[0][1,None,None]])
# Full slice
X = DummyNode([np.random.randn(3,4)],
plates=(3,4),
dims=((),))
Y = X[:,2]
self.assertMessageToChild(Y, [X.u[0][:,2]])
# Slice with start
X = DummyNode([np.random.randn(3,4)],
plates=(3,4),
dims=((),))
Y = X[1:,2]
self.assertMessageToChild(Y, [X.u[0][1:,2]])
# Slice with end
X = DummyNode([np.random.randn(3,4)],
plates=(3,4),
dims=((),))
Y = X[:2,2]
self.assertMessageToChild(Y, [X.u[0][:2,2]])
# Slice with step
X = DummyNode([np.random.randn(3,4)],
plates=(3,4),
dims=((),))
Y = X[::2,2]
self.assertMessageToChild(Y, [X.u[0][::2,2]])
# Slice for broadcasted axes
X = DummyNode([np.random.randn(4)],
plates=(3,4),
dims=((),))
Y = X[0:2:2,2]
self.assertMessageToChild(Y, [X.u[0][2]])
X = DummyNode([np.random.randn(1,4)],
plates=(3,4),
dims=((),))
Y = X[0:2:2,2]
self.assertMessageToChild(Y, [X.u[0][0:1,2]])
# One ellipsis
X = DummyNode([np.random.randn(3,4)],
plates=(3,4),
dims=((),))
Y = X[...,2]
self.assertMessageToChild(Y, [X.u[0][...,2]])
# Ellipsis over broadcasted axes
X = DummyNode([np.random.randn(5,6)],
plates=(3,4,5,6),
dims=((),))
Y = X[1,...,2]
self.assertMessageToChild(Y, [X.u[0][:,2]])
X = DummyNode([np.random.randn(3,1,5,6)],
plates=(3,4,5,6),
dims=((),))
Y = X[1,...,2]
self.assertMessageToChild(Y, [X.u[0][1,:,:,2]])
# Multiple ellipsis
X = DummyNode([np.random.randn(2,3,4,5)],
plates=(2,3,4,5),
dims=((),))
Y = X[...,2,...]
self.assertMessageToChild(Y, [X.u[0][:,:,2,:]])
# Ellipsis when dimensions
X = DummyNode([np.random.randn(2,3,4)],
plates=(2,3),
dims=((4,),))
Y = X[...,2]
self.assertMessageToChild(Y, [X.u[0][:,2,:]])
# Indexing for multiple moments
X = DummyNode([np.random.randn(2,3,4),
np.random.randn(2,3)],
plates=(2,3),
dims=((4,),()))
Y = X[1,1]
self.assertMessageToChild(Y, [X.u[0][1,1],
X.u[1][1,1]])
pass
0
Example 95
def test_message_to_parent(self):
"""
Test message to parent of X[..] node operator.
"""
class ParentNode(Node):
_moments = Moments()
_parent_moments = ()
def _get_id_list(self):
return []
class ChildNode(Node):
_moments = Moments()
_parent_moments = (Moments(),)
def __init__(self, X, m, mask, **kwargs):
super().__init__(X, **kwargs)
self.m = m
self.mask2 = mask
def _message_to_parent(self, index, u_parent=None):
return self.m
def _mask_to_parent(self, index):
return self.mask2
def _get_id_list(self):
return []
# General broadcasting
V = ParentNode(plates=(3,3,3),
dims=((),))
X = V[...]
m = [ np.random.randn(3,1) ]
msg = [ np.zeros((1,3,1)) ]
msg[0][:,:,:] = m[0]
Y = ChildNode(X, m, True, dims=((),))
X._update_mask()
self.assertMessage(X._message_to_parent(0),
msg)
# Integer indices
V = ParentNode(plates=(3,4),
dims=((),))
X = V[2,1]
m = [np.random.randn()]
msg = [ np.zeros((3,4)) ]
msg[0][2,1] = m[0]
Y = ChildNode(X, m, True, dims=((),))
X._update_mask()
self.assertMessage(X._message_to_parent(0),
msg)
# Integer indices with broadcasting
V = ParentNode(plates=(3,3),
dims=((),))
X = V[2,2]
m = [ np.random.randn(1) ]
msg = [ np.zeros((3,3)) ]
msg[0][2,2] = m[0]
Y = ChildNode(X, m, True, dims=((),))
X._update_mask()
self.assertMessage(X._message_to_parent(0),
msg)
# Slice indices
V = ParentNode(plates=(2,3,4,5),
dims=((),))
X = V[:,:2,1:,::2]
m = [np.random.randn(2,2,3,3)]
msg = [ np.zeros((2,3,4,5)) ]
msg[0][:,:2,1:,::2] = m[0]
Y = ChildNode(X, m, True, dims=((),))
X._update_mask()
self.assertMessage(X._message_to_parent(0),
msg)
# Full slice with broadcasting
V = ParentNode(plates=(2,3),
dims=((),))
X = V[:,:]
m = [np.random.randn(1)]
msg = [ np.zeros((1,1)) ]
msg[0][:] = m[0]
Y = ChildNode(X, m, True, dims=((),))
X._update_mask()
self.assertMessage(X._message_to_parent(0),
msg)
# Start slice with broadcasting
V = ParentNode(plates=(3,3,3,3),
dims=((),))
X = V[0:,1:,-2:,-3:]
m = [np.random.randn(1,1)]
msg = [ np.zeros((1,3,3,1)) ]
msg[0][:,1:,-2:,:] = m[0]
Y = ChildNode(X, m, True, dims=((),))
X._update_mask()
self.assertMessage(X._message_to_parent(0),
msg)
# End slice with broadcasting
V = ParentNode(plates=(3,3,3,3),
dims=((),))
X = V[:2,:3,:4,:-1]
m = [np.random.randn(1,1)]
msg = [ np.zeros((3,1,1,3)) ]
msg[0][:2,:,:,:-1] = m[0]
Y = ChildNode(X, m, True, dims=((),))
X._update_mask()
self.assertMessage(X._message_to_parent(0),
msg)
# Step slice with broadcasting
V = ParentNode(plates=(3,3,1),
dims=((),))
X = V[::1,::2,::2]
m = [np.random.randn(1)]
msg = [ np.zeros((1,3,1)) ]
msg[0][:,::2,:] = m[0]
Y = ChildNode(X, m, True, dims=((),))
X._update_mask()
self.assertMessage(X._message_to_parent(0),
msg)
# Ellipsis
V = ParentNode(plates=(3,3,3),
dims=((),))
X = V[...,0]
m = [np.random.randn(3,3)]
msg = [ np.zeros((3,3,3)) ]
msg[0][:,:,0] = m[0]
Y = ChildNode(X, m, True, dims=((),))
X._update_mask()
self.assertMessage(X._message_to_parent(0),
msg)
# New axes
V = ParentNode(plates=(3,3),
dims=((),))
X = V[None,:,None,None,:]
m = [np.random.randn(1,3,1,1,3)]
msg = [ np.zeros((3,3)) ]
msg[0][:,:] = m[0][0,:,0,0,:]
Y = ChildNode(X, m, True, dims=((),))
X._update_mask()
self.assertMessage(X._message_to_parent(0),
msg)
# New axes with broadcasting
V = ParentNode(plates=(3,3),
dims=((),))
X = V[None,:,None,:,None]
m = [np.random.randn(1,3,1)]
msg = [ np.zeros((1,3)) ]
msg[0][:,:] = m[0][0,:,0]
Y = ChildNode(X, m, True, dims=((),))
X._update_mask()
self.assertMessage(X._message_to_parent(0),
msg)
# Multiple messages
V = ParentNode(plates=(3,),
dims=((),()))
X = V[:]
m = [np.random.randn(3),
np.random.randn(3)]
msg = [ np.zeros((3)), np.zeros((3)) ]
msg[0][:] = m[0][:]
msg[1][:] = m[1][:]
Y = ChildNode(X, m, True, dims=((),()))
X._update_mask()
self.assertMessage(X._message_to_parent(0),
msg)
# Non-scalar variables
V = ParentNode(plates=(2,3),
dims=((4,),))
X = V[...]
m = [np.random.randn(2,3,4)]
msg = [ np.zeros((2,3,4)) ]
msg[0][:,:,:] = m[0][:,:,:]
Y = ChildNode(X, m, True, dims=((4,),))
X._update_mask()
self.assertMessage(X._message_to_parent(0),
msg)
# Missing values
V = ParentNode(plates=(3,3,3),
dims=((3,),))
X = V[:,0,::2,None]
m = [np.random.randn(3,2,1,3)]
# mask shape: (3, 2, 1)
mask = np.array([ [[True], [False]],
[[False], [False]],
[[False], [True]] ])
msg = [ np.zeros((3,3,3,3)) ]
msg[0][:,0,::2,:] = (m[0] * mask[...,None])[:,:,0,:]
Y = ChildNode(X, m, mask, dims=((3,),))
X._update_mask()
self.assertMessage(X._message_to_parent(0),
msg)
# Found bug: int index after slice
V = ParentNode(plates=(3,3),
dims=((),))
X = V[:,0]
m = [np.random.randn(3)]
msg = [ np.zeros((3,3)) ]
msg[0][:,0] = m[0]
Y = ChildNode(X, m, True, dims=((),))
X._update_mask()
self.assertMessage(X._message_to_parent(0),
msg)
# Found bug: message requires reshaping after reverse indexing
pass
0
Example 96
Project: warpedLMM Source File: linalg.py
def ppca_missing_data_at_random(Y, Q, iters=100):
"""
EM implementation of Probabilistic pca for when there is missing data.
Taken from <SheffieldML, https://github.com/SheffieldML>
.. math:
\\mathbf{Y} = \mathbf{XW} + \\epsilon \\text{, where}
\\epsilon = \\mathcal{N}(0, \\sigma^2 \mathbf{I})
:returns: X, W, sigma^2
"""
from numpy.ma import dot as madot
import diag
from GPy.util.subarray_and_sorting import common_subarrays
import time
debug = 1
# Initialise W randomly
N, D = Y.shape
W = np.random.randn(Q, D) * 1e-3
Y = np.ma.masked_invalid(Y, copy=1)
nu = 1.
#num_obs_i = 1./Y.count()
Ycentered = Y - Y.mean(0)
X = np.zeros((N,Q))
cs = common_subarrays(Y.mask)
cr = common_subarrays(Y.mask, 1)
Sigma = np.zeros((N, Q, Q))
Sigma2 = np.zeros((N, Q, Q))
mu = np.zeros(D)
"""
if debug:
import matplotlib.pyplot as pylab
fig = pylab.figure("FIT MISSING DATA");
ax = fig.gca()
ax.cla()
lines = pylab.plot(np.zeros((N,Q)).dot(W))
"""
W2 = np.zeros((Q,D))
for i in range(iters):
# Sigma = np.linalg.solve(diag.add(madot(W,W.T), nu), diag.times(np.eye(Q),nu))
# exp_x = madot(madot(Ycentered, W.T),Sigma)/nu
# Ycentered = (Y - exp_x.dot(W).mean(0))
# #import ipdb;ipdb.set_trace()
# #Ycentered = mu
# W = np.linalg.solve(madot(exp_x.T,exp_x) + Sigma, madot(exp_x.T, Ycentered))
# nu = (((Ycentered - madot(exp_x, W))**2).sum(0) + madot(W.T,madot(Sigma,W)).sum(0)).sum()/N
for csi, (mask, index) in enumerate(cs.iteritems()):
mask = ~np.array(mask)
Sigma2[index, :, :] = nu * np.linalg.inv(diag.add(W2[:,mask].dot(W2[:,mask].T), nu))
#X[index,:] = madot((Sigma[csi]/nu),madot(W,Ycentered[index].T))[:,0]
X2 = ((Sigma2/nu) * (madot(Ycentered,W2.T).base)[:,:,None]).sum(-1)
mu2 = (Y - X.dot(W)).mean(0)
for n in range(N):
Sigma[n] = nu * np.linalg.inv(diag.add(W[:,~Y.mask[n]].dot(W[:,~Y.mask[n]].T), nu))
X[n, :] = (Sigma[n]/nu).dot(W[:,~Y.mask[n]].dot(Ycentered[n,~Y.mask[n]].T))
for d in range(D):
mu[d] = (Y[~Y.mask[:,d], d] - X[~Y.mask[:,d]].dot(W[:, d])).mean()
Ycentered = (Y - mu)
nu3 = 0.
for cri, (mask, index) in enumerate(cr.iteritems()):
mask = ~np.array(mask)
W2[:,index] = np.linalg.solve(X[mask].T.dot(X[mask]) + Sigma[mask].sum(0), madot(X[mask].T, Ycentered[mask,index]))[:,None]
W2[:,index] = np.linalg.solve(X.T.dot(X) + Sigma.sum(0), madot(X.T, Ycentered[:,index]))
#nu += (((Ycentered[mask,index] - X[mask].dot(W[:,index]))**2).sum(0) + W[:,index].T.dot(Sigma[mask].sum(0).dot(W[:,index])).sum(0)).sum()
nu3 += (((Ycentered[index] - X.dot(W[:,index]))**2).sum(0) + W[:,index].T.dot(Sigma.sum(0).dot(W[:,index])).sum(0)).sum()
nu3 /= N
nu = 0.
nu2 = 0.
W = np.zeros((Q,D))
for j in range(D):
W[:,j] = np.linalg.solve(X[~Y.mask[:,j]].T.dot(X[~Y.mask[:,j]]) + Sigma[~Y.mask[:,j]].sum(0), madot(X[~Y.mask[:,j]].T, Ycentered[~Y.mask[:,j],j]))
nu2f = np.tensordot(W[:,j].T, Sigma[~Y.mask[:,j],:,:], [0,1]).dot(W[:,j])
nu2s = W[:,j].T.dot(Sigma[~Y.mask[:,j],:,:].sum(0).dot(W[:,j]))
nu2 += (((Ycentered[~Y.mask[:,j],j] - X[~Y.mask[:,j],:].dot(W[:,j]))**2) + nu2f).sum()
for i in range(N):
if not Y.mask[i,j]:
nu += ((Ycentered[i,j] - X[i,:].dot(W[:,j]))**2) + W[:,j].T.dot(Sigma[i,:,:].dot(W[:,j]))
nu /= N
nu2 /= N
nu4 = (((Ycentered - X.dot(W))**2).sum(0) + W.T.dot(Sigma.sum(0).dot(W)).sum(0)).sum()/N
import ipdb;ipdb.set_trace()
"""
if debug:
#print Sigma[0]
print "nu:", nu, "sum(X):", X.sum()
pred_y = X.dot(W)
for x, l in zip(pred_y.T, lines):
l.set_ydata(x)
ax.autoscale_view()
ax.set_ylim(pred_y.min(), pred_y.max())
fig.canvas.draw()
time.sleep(.3)
"""
return np.asarray_chkfinite(X), np.asarray_chkfinite(W), nu
0
Example 97
def __init__(self, n=[301, 501], select=["all", ], n_max_bs=700,
n_max_slow=700, transform_repeat=1):
"""
Benchmark performance of different iAbel/fAbel implementations.
Parameters
----------
n: integer
a list of arrays sizes for the benchmark
(assuming 2D square arrays (n,n))
select: list of str
list of transforms to benchmark select=['all',] (default) or
choose transforms:
select=['basex', 'direct_Python', 'direct_C', 'hansenlaw',
'linbasex' 'onion_bordas, 'onion_peeling', 'two_point',
'three_point']
n_max_bs: integer
since the basis sets generation takes a long time,
do not run this benchmark for implementations that use basis sets
for n > n_max_bs
n_max_slow: integer
maximum n run for the "slow" transform methods, so far including
only the "direct_python" implementation.
"""
from timeit import Timer
import time
self.n = n
transform = {
'basex': basex.basex_core_transform,
'basex_bs': basex.get_bs_basex_cached,
'direct_Python': direct.direct_transform,
'direct_C': direct.direct_transform,
'hansenlaw': hansenlaw.hansenlaw_transform,
'linbasex': linbasex._linbasex_transform_with_basis,
'linbasex_bs': linbasex._bs_linbasex,
'onion_bordas': onion_bordas.onion_bordas_transform,
'onion_peeling': dasch.dasch_transform,
'onion_peeling_bs': dasch._bs_onion_peeling,
'two_point': dasch.dasch_transform,
'two_point_bs': dasch._bs_two_point,
'three_point': dasch.dasch_transform,
'three_point_bs': dasch._bs_three_point,
}
# result dicts
res = {}
res['bs'] = {'basex_bs': [], 'linbasex_bs': [], 'onion_peeling_bs': [],
'two_point_bs': [], 'three_point_bs': []}
res['forward'] = {'direct_Python': [], 'hansenlaw': []}
res['inverse'] = {'basex': [], 'direct_Python': [], 'hansenlaw': [],
'linbasex': [],
'onion_bordas': [], 'onion_peeling': [],
'two_point': [], 'three_point': []}
if direct.cython_ext:
res['forward']['direct_C'] = []
res['inverse']['direct_C'] = []
# delete all keys not present in 'select' input parameter
if "all" not in select:
for trans in select:
if trans not in res['inverse'].keys():
raise ValueError("'{}' is not a valid transform method".
format(trans), res['inverse'].keys())
for direction in ['forward', 'inverse']:
rm = []
for abel in res[direction]:
if abel not in select:
rm.append(abel)
for x in rm:
del res[direction][x]
# repeat for 'bs' which has append '_bs'
rm = []
for abel in res['bs']:
if abel[:-3] not in select:
rm.append(abel)
for x in rm:
del res['bs'][x]
# ---- timing tests for various image sizes nxn
for ni in n:
ni = int(ni)
x = np.random.randn(ni, ni)
# basis set evaluation --------------
basis = {}
for method in res['bs'].keys():
if method[:-3] == 'basex': # special case
if ni <= n_max_bs:
# calculate and store basex basis matrix
t = time.time()
basis[method[:-3]] = transform[method](ni, ni,
basis_dir=None)
res['bs'][method].append((time.time()-t)*1000)
else:
basis[method[:-3]] = None,
res['bs'][method].append(np.nan)
else:
# calculate and store basis matrix
t = time.time()
# store basis calculation. NB a tuple to accomodate basex
basis[method[:-3]] = transform[method](ni),
res['bs'][method].append((time.time()-t)*1000)
# Abel transforms ---------------
for cal in ["forward", "inverse"]:
for method in res[cal].keys():
if method in basis.keys():
if basis[method][0] is not None:
# have basis calculation
res[cal][method].append(Timer(
lambda: transform[method](x, *basis[method])).
timeit(number=transform_repeat)*1000/
transform_repeat)
else:
# no calculation available
res[cal][method].append(np.nan)
elif method[:6] == 'direct': # special case 'direct'
if method[7] == 'P' and (ni > n_max_slow):
res[cal][method].append(np.nan)
else:
res[cal][method].append(Timer(
lambda: transform[method](x, backend=method[7:],
direction=cal)).
timeit(number=transform_repeat)*1000/
transform_repeat)
else:
# full calculation for everything else
res[cal][method].append(Timer(
lambda: transform[method](x, direction=cal)).
timeit(number=transform_repeat)*1000/
transform_repeat)
self.fabel = res['forward']
self.bs = res['bs']
self.iabel = res['inverse']
0
Example 98
Project: pyFFTW Source File: test_pyfftw_call.py
def test_call_with_unaligned(self):
'''Make sure the right thing happens with unaligned data.
'''
input_array = (numpy.random.randn(*self.input_array.shape)
+ 1j*numpy.random.randn(*self.input_array.shape))
output_array = self.fft(
input_array=byte_align(input_array.copy(), n=16)).copy()
input_array = byte_align(input_array, n=16)
output_array = byte_align(output_array, n=16)
# Offset by one from 16 byte aligned to guarantee it's not
# 16 byte aligned
a = byte_align(input_array.copy(), n=16)
a__ = empty_aligned(numpy.prod(a.shape)*a.itemsize+1, dtype='int8',
n=16)
a_ = a__[1:].view(dtype=a.dtype).reshape(*a.shape)
a_[:] = a
# Create a different second array the same way
b = byte_align(output_array.copy(), n=16)
b__ = empty_aligned(numpy.prod(b.shape)*a.itemsize+1, dtype='int8',
n=16)
b_ = b__[1:].view(dtype=b.dtype).reshape(*b.shape)
b_[:] = a
# Set up for the first array
fft = FFTW(input_array, output_array)
a_[:] = a
output_array = fft().copy()
# Check a_ is not aligned...
self.assertRaisesRegex(ValueError, 'Invalid input alignment',
self.fft.update_arrays, *(a_, output_array))
# and b_ too
self.assertRaisesRegex(ValueError, 'Invalid output alignment',
self.fft.update_arrays, *(input_array, b_))
# But it should still work with the a_
fft(a_)
# However, trying to update the output will raise an error
self.assertRaisesRegex(ValueError, 'Invalid output alignment',
self.fft.update_arrays, *(input_array, b_))
# Same with SIMD off
fft = FFTW(input_array, output_array, flags=('FFTW_UNALIGNED',))
fft(a_)
self.assertRaisesRegex(ValueError, 'Invalid output alignment',
self.fft.update_arrays, *(input_array, b_))
0
Example 99
Project: attention-lvcsr Source File: test_gradient.py
def test_subgraph_grad():
# Tests that the grad method with no known_grads
# matches what happens if you use successive subgraph_grads
x = theano.tensor.fvector('x')
t = theano.tensor.fvector('t')
w1 = theano.shared(np.random.randn(3, 4))
w2 = theano.shared(np.random.randn(4, 2))
a1 = theano.tensor.tanh(theano.tensor.dot(x, w1))
a2 = theano.tensor.tanh(theano.tensor.dot(a1, w2))
cost2 = theano.tensor.sqr(a2 - t).sum()
cost2 += theano.tensor.sqr(w2.sum())
cost1 = theano.tensor.sqr(w1.sum())
params = [[w2], [w1]]
costs = [cost2, cost1]
grad_ends = [[a1], [x]]
inputs = [t, x]
rng = np.random.RandomState([2012, 11, 15])
values = [rng.randn(2), rng.randn(3)]
values = [np.cast[ipt.dtype](value) for ipt, value in zip(inputs, values)]
wrt = [w2, w1]
cost = cost2 + cost1
true_grads = theano.grad(cost, wrt)
true_grads = theano.function(inputs, true_grads)
true_grads = true_grads(*values)
next_grad = None
param_grads = []
for i in xrange(2):
param_grad, next_grad = theano.subgraph_grad(
wrt=params[i], end=grad_ends[i],
start=next_grad, cost=costs[i]
)
next_grad = OrderedDict(izip(grad_ends[i], next_grad))
param_grads.extend(param_grad)
pgrads = theano.function(inputs, param_grads)
pgrads = pgrads(*values)
for true_grad, pgrad in zip(true_grads, pgrads):
assert(np.sum(np.abs(true_grad - pgrad)) < 0.00001)
0
Example 100
def setUp(self):
"""
Setup the dataframes used for the groupby tests derived from pandas
"""
self.date_rng = bdate_range('1/1/2005', periods=250)
self.string_idx = Index([rands(8).upper() for x in range(250)])
self.group_id = Series([x[0] for x in self.string_idx],
index=self.string_idx)
self.group_dict = dict((key, value) for key, value in
compat.iteritems(self.group_id))
self.col_idx = Index(['A', 'B', 'C', 'D', 'E'])
rand_matrix = np.random.randn(250, 5)
self.string_matrix = DataFrame(rand_matrix, columns=self.col_idx,
index=self.string_idx)
self.time_matrix = DataFrame(rand_matrix, columns=self.col_idx,
index=self.date_rng)
self.time_series = tm.makeTimeSeries()
self.seriesd = tm.getSeriesData()
self.tsd = tm.getTimeSeriesData()
self.frame = DataFrame(self.seriesd)
self.tsframe = DataFrame(self.tsd)
self.pd_df_foobar = DataFrame({'A': ['foo', 'bar', 'foo', 'bar',
'foo', 'bar', 'foo', 'foo'],
'B': ['one', 'one', 'two', 'three',
'two', 'two', 'one', 'three'],
'C': np.random.randn(8),
'D': np.random.randn(8)})
self.df_mixed_floats = DataFrame({'A': ['foo', 'bar', 'foo', 'bar',
'foo', 'bar', 'foo', 'foo'],
'B': ['one', 'one', 'two', 'three',
'two', 'two', 'one', 'three'],
'C': np.random.randn(8),
'D': np.array(np.random.randn(8),
dtype='float32')})
index = MultiIndex(levels=[['foo', 'bar', 'baz', 'qux'],
['one', 'two', 'three']],
labels=[[0, 0, 0, 1, 1, 2, 2, 3, 3, 3],
[0, 1, 2, 0, 1, 1, 2, 0, 1, 2]],
names=['first', 'second'])
self.mframe = DataFrame(np.random.randn(10, 3), index=index,
columns=['A', 'B', 'C'])
self.three_group = DataFrame({'A': ['foo', 'foo', 'foo', 'foo',
'bar', 'bar', 'bar', 'bar',
'foo', 'foo', 'foo'],
'B': ['one', 'one', 'one', 'two',
'one', 'one', 'one', 'two',
'two', 'two', 'one'],
'C': ['dull', 'dull', 'shiny', 'dull',
'dull', 'shiny', 'shiny', 'dull',
'shiny', 'shiny', 'shiny'],
'D': np.random.randn(11),
'E': np.random.randn(11),
'F': np.random.randn(11)})
super(self.__class__, self).setUp()