numpy.random.randn

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 7

Example 51

Project: pylearn2 Source File: autoencoder.py
Function: init
    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]

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'

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))

Example 54

Project: mlp-classifier Source File: neuralnet.py
Function: initialize_weights
	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:])]

Example 55

Project: deepgaze Source File: motion_tracking.py
Function: predict
    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

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)

Example 57

Project: nideep Source File: test_balance.py
Function: set_up
    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

Example 58

Project: matplotlib2tikz Source File: histogram.py
Function: plot
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

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

Example 60

Project: downhill Source File: dataset_test.py
Function: test_batch_size
    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)

Example 61

Project: lightning-python Source File: test_client.py
Function: test_create
    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')

Example 62

Project: nupic.research Source File: sparse_net.py
Function: reset
  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 = {}

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)

Example 64

Project: dl4mt-tutorial Source File: lm.py
Function: norm_weight
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')

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__)

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

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]

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)

Example 69

Project: spectralDNS Source File: test_shentransforms.py
Function: test_helmholtz
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

Example 70

Project: filterpy Source File: test_sqrtkf.py
Function: test_noisy_1d
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()

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

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')

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

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()

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()

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()

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)

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()

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()

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

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

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)

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))

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)

Example 85

Project: simpeg Source File: Tests.py
Function: check_derivative
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

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)

Example 87

Project: QuantEcon.py Source File: lss.py
Function: simulate
    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

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

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

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)

Example 91

Project: bayespy Source File: test_gaussian_markov_chain.py
Function: test_message_to_child
    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

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

Example 93

Project: cortex Source File: tsne.py
Function: tsne
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;

Example 94

Project: bayespy Source File: test_node.py
Function: test_message_to_child
    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

Example 95

Project: bayespy Source File: test_node.py
Function: test_message_to_parent
    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

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

Example 97

Project: PyAbel Source File: benchmark.py
Function: init
    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']

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_))

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)

Example 100

Project: sparklingpandas Source File: pandas_groupby_tests.py
Function: set_up
    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()
See More Examples - Go to Next Page
Page 1 Page 2 Selected Page 3 Page 4