def _train(self, x, y):
        """
        :param x: Array of different input observations.
        :type x: numpy.ndarray

        :param y: Array of size (x.shape[0], output_dim) that contains the 
            observed output to the input x's.
        :type y: numpy.ndarray
        """
        # initialize internal vars if necessary
        if self._xTx is None:
            if self.with_bias:
                x_size = self._input_dim + 1
            else:
                x_size = self._input_dim
            self._xTx = numx.zeros((x_size, x_size), self._dtype)
            self._xTy = numx.zeros((x_size, self._output_dim), self._dtype)

        if self.with_bias:
            x = self._add_constant(x)

        # update internal variables
        self._xTx += mult(x.T, x)
        self._xTy += mult(x.T, y)
        self._tlen += x.shape[0]

 def _init_internals(self):
     input_dim = self.input_dim
     self._mean = numx.zeros((input_dim,), dtype='d')
     self._var = numx.zeros((input_dim,), dtype='d')
     self._tlen = 0
     self._diff2 = numx.zeros((input_dim,), dtype='d')
     self._initialized = 1

    def _execute(self, x):
        degree = self._degree
        dim = self.input_dim
        n = x.shape[1]

        # preallocate memory
        dexp = numx.zeros((self.output_dim, x.shape[0]), dtype=self.dtype)
        # copy monomials of degree 1
        dexp[0:n, :] = x.T

        k = n
        prec_end = 0
        next_lens = numx.ones((dim+1, ))
        next_lens[0] = 0
        for i in range(2, degree+1):
            prec_start = prec_end
            prec_end += nmonomials(i-1, dim)
            prec = dexp[prec_start:prec_end, :]

            lens = next_lens[:-1].cumsum(axis=0)
            next_lens = numx.zeros((dim+1, ))
            for j in range(dim):
                factor = prec[lens[j]:, :]
                len_ = factor.shape[0]
                dexp[k:k+len_, :] = x[:, j] * factor
                next_lens[j+1] = len_
                k = k+len_

        return dexp.T

 def __init__(self, n, d, real_dtype="d", integer_dtype="l"):
     """Initializes an object of type 'OneDimensionalHitParade'.
     
     :param n: Number of maxima and minima to remember.
     :type n: int
     
     :param d: Minimum gap between two hits.
     :type d: int
     
     :param real_dtype: Datatype of sequence items
     :type real_dtype: numpy.dtype or str
     
     :param integer_dtype: Datatype of sequence indices
     :type integer_dtype: numpy.dtype or str
     """
     self.n = int(n)
     self.d = int(d)
     self.iM = numx.zeros((n, ), dtype=integer_dtype)
     self.im = numx.zeros((n, ), dtype=integer_dtype)
     
     real_dtype = numx.dtype(real_dtype)
     if real_dtype in mdp.utils.get_dtypes('AllInteger'):
         max_num = numx.iinfo(real_dtype).max
         min_num = numx.iinfo(real_dtype).min
     else:
         max_num = numx.finfo(real_dtype).max
         min_num = numx.finfo(real_dtype).min
     self.M = numx.array([min_num]*n, dtype=real_dtype)
     self.m = numx.array([max_num]*n, dtype=real_dtype)
     
     self.lM = 0
     self.lm = 0

 def test_mixed_dict(self):
     """Test msg being a dict containing an array."""
     rescont = MessageResultContainer()
     msg1 = {
         "f": 2,
         "a": np.zeros((10,3), 'int'),
         "b": "aaa",
         "c": 1,
     }
     msg2 = {
         "a": np.ones((15,3), 'int'),
         "b": "bbb",
         "c": 3,
         "d": 1,
     }
     rescont.add_message(msg1)
     rescont.add_message(msg2)
     combined_msg = rescont.get_message()
     a = np.zeros((25,3), 'int')
     a[10:] = 1
     reference_msg = {"a": a, "c": 4, "b": "aaabbb", "d": 1, "f": 2}
     assert np.all(reference_msg["a"] == reference_msg["a"])
     combined_msg.pop("a")
     reference_msg.pop("a")
     assert combined_msg == reference_msg

    def __init__(self, n, d, real_dtype="d", integer_dtype="l"):
        """
        Input arguments:
        n -- Number of maxima and minima to remember
        d -- Minimum gap between two hits

        real_dtype -- dtype of sequence items
        integer_dtype -- dtype of sequence indices
        Note: be careful with dtypes!
        """
        self.n = int(n)
        self.d = int(d)
        self.iM = numx.zeros((n, ), dtype=integer_dtype)
        self.im = numx.zeros((n, ), dtype=integer_dtype)
        
        real_dtype = numx.dtype(real_dtype)
        if real_dtype in mdp.utils.get_dtypes('AllInteger'):
            max_num = numx.iinfo(real_dtype).max
            min_num = numx.iinfo(real_dtype).min
        else:
            max_num = numx.finfo(real_dtype).max
            min_num = numx.finfo(real_dtype).min
        self.M = numx.array([min_num]*n, dtype=real_dtype)
        self.m = numx.array([max_num]*n, dtype=real_dtype)
        
        self.lM = 0
        self.lm = 0

 def test_incompatible_arrays(self):
     """Test with incompatible arrays."""
     rescont = MessageResultContainer()
     msgs = [{"a":  np.zeros((10,3))}, {"a":  np.zeros((10,4))}]
     for msg in msgs:
         rescont.add_message(msg)
     pytest.raises(ValueError, rescont.get_message)

 def output_sizes(self, n):
     """Return the individual output sizes of each expansion function
     when the input has lenght n."""
     sizes = numx.zeros(len(self.funcs))
     x = numx.zeros((1,n))
     for i, func in enumerate(self.funcs):
         outx = func(x)
         sizes[i] = outx.shape[1]
     return sizes

 def _get_rnd_permutation(self, dim):
     # return a random permut matrix with the right dimensions and type
     zero = numx.zeros((dim, dim), dtype=self.dtype)
     row = numx_rand.permutation(dim)
     for col in range(dim):
         zero[row[col], col] = 1.
     return zero

 def output_sizes(self, n):
     """Return the individual output sizes of each expansion function
     when the input has lenght n.
     
     :param n: Input dimension,
     :type: int
     
     :return: The individual output sizes of each expansion function.
     :rtype: list
     """
     sizes = numx.zeros(len(self.funcs), dtype=numx.int64)
     x = numx.zeros((1,n))
     for i, func in enumerate(self.funcs):
         outx = func(x)
         sizes[i] = outx.shape[1]
     return sizes

    def _execute(self, x):
        assert x.shape[0] == 1

        if self.sliding_wnd == None:
            self._init_sliding_window()

        gap = self.gap
        rows = self.sliding_wnd.shape[0]
        cols = self.output_dim
        n = self.input_dim

        new_row = numx.zeros(cols, dtype=self.dtype)
        new_row[:n] = x

        # Slide
        if self.slide:
            self.sliding_wnd[:-1, :] = self.sliding_wnd[1:, :]

        # Delay
        if self.cur_idx-gap >= 0:
            new_row[n:] = self.sliding_wnd[self.cur_idx-gap, :-n]

        # Add new row to matrix
        self.sliding_wnd[self.cur_idx, :] = new_row

        if self.cur_idx < rows-1:
            self.cur_idx = self.cur_idx+1 
        else:
            self.slide = True

        return new_row[numx.newaxis,:]

    def get_quadratic_form(self, nr):
        """
        Return the matrix H, the vector f and the constant c of the
        quadratic form 1/2 x'Hx + f'x + c that defines the output
        of the component 'nr' of the SFA node.
        """
        if self.sf is None:
            self._if_training_stop_training()

        sf = self.sf[:, nr]
        c = -mult(self.avg, sf)
        n = self.input_dim
        f = sf[:n]
        h = numx.zeros((n, n), dtype=self.dtype)
        k = n
        for i in range(n):
            for j in range(n):
                if j > i:
                    h[i, j] = sf[k]
                    k = k+1
                elif j == i:
                    h[i, j] = 2*sf[k]
                    k = k+1
                else:
                    h[i, j] = h[j, i]

        return QuadraticForm(h, f, c, dtype=self.dtype)

    def _sample_v(self, h, sample_l=False, concatenate=True):
        # returns  P(v=1|h,W,b), a sample from it, P(l=1|h,W,b),
        # and a sample from it

        ldim, vdim = self._labels_dim, self._visible_dim

        # activation
        a = self.bv + mult(h, self.w.T)
        av, al = a[:, :vdim], a[:, vdim:]

        # ## visible units: logistic activation
        probs_v = old_div(1.,(1. + exp(-av)))
        v = (probs_v > random(probs_v.shape)).astype('d')

        # ## label units: softmax activation
        # subtract maximum to regularize exponent
        exponent = al - rrep(al.max(axis=1), ldim)
        probs_l = exp(exponent)
        probs_l /= rrep(probs_l.sum(axis=1), ldim)

        if sample_l:
            # ?? todo: I'm sure this can be optimized
            l = numx.zeros((h.shape[0], ldim))
            for t in range(h.shape[0]):
                l[t, :] = mdp.numx_rand.multinomial(1, probs_l[t, :])
        else:
            l = probs_l.copy()

        if concatenate:
            probs = numx.concatenate((probs_v, probs_l), axis=1)
            x = numx.concatenate((v, l), axis=1)
            return probs, x
        else:
            return probs_v, probs_l, v, l

    def _execute(self, x):
        #----------------------------------------------------
        # similar algorithm to that within self.stop_training()
        #  refer there for notes & comments on code
        #----------------------------------------------------
        N = self.data.shape[0]
        Nx = x.shape[0]
        W = numx.zeros((Nx, N), dtype=self.dtype)

        k, r = self.k, self.r
        d_out = self.output_dim
        Q_diag_idx = numx.arange(k)

        for row in range(Nx):
            #find nearest neighbors of x in M
            M_xi = self.data-x[row]
            nbrs = numx.argsort( (M_xi**2).sum(1) )[:k]
            M_xi = M_xi[nbrs]

            #find corrected covariance matrix Q
            Q = mult(M_xi, M_xi.T)
            if r is None and k > d_out:
                sig2 = (svd(M_xi, compute_uv=0))**2
                r = numx.sum(sig2[d_out:])
                Q[Q_diag_idx, Q_diag_idx] += r
            if r is not None:
                Q[Q_diag_idx, Q_diag_idx] += r

            #solve for weights
            w = self._refcast(numx_linalg.solve(Q , numx.ones(k)))
            w /= w.sum()
            W[row, nbrs] = w

        #multiply weights by result of SVD from training
        return numx.dot(W, self.training_projection)

 def _update_mean(self, x, label):
     """Update the mean with data for a single label."""
     if label not in self.label_means:
         self.label_means[label] = numx.zeros(self.input_dim)
         self.n_label_samples[label] = 0
     # TODO: use smarter summing to avoid rounding errors
     self.label_means[label] += numx.sum(x, axis=0)
     self.n_label_samples[label] += len(x)

 def _get_contrast(self, covs, bica_bsfa = None):
     if bica_bsfa is None:
         bica_bsfa = self._bica_bsfa
     # return current value of the contrast
     R = self.output_dim
     ncovs = covs.ncovs
     covs = covs.covs
     icaweights = self.icaweights
     sfaweights = self.sfaweights
     # unpack the bsfa and bica coefficients
     bica, bsfa = bica_bsfa
     sfa = numx.zeros((ncovs, ), dtype=self.dtype)
     ica = numx.zeros((ncovs, ), dtype=self.dtype)
     for t in range(ncovs):
         sq_corr =  covs[:R, :R, t]*covs[:R, :R, t]
         sfa[t] = sq_corr.trace()
         ica[t] = 2*_triu(sq_corr, 1).ravel().sum()
     return (bsfa*sfaweights*sfa).sum(), (bica*icaweights*ica).sum()

 def test_execute_routing(self):
     """Test the standard routing for messages."""
     sboard = BiSwitchboard(input_dim=3, connections=[2,0,1])
     x = n.array([[1,2,3],[4,5,6]])
     msg = {
         "string": "blabla",
         "list": [1,2],
         "data": x.copy(),  # should be mapped by switchboard
         "data2": n.zeros(3),  # should not be modified
         "data3": n.zeros((3,4)),  # should not be modified
     }
     y, out_msg = sboard.execute(x, msg)
     reference_y = n.array([[3,1,2],[6,4,5]])
     assert (y == reference_y).all()
     assert out_msg["string"] == msg["string"]
     assert out_msg["list"] == msg["list"]
     assert n.all(out_msg["data"] == reference_y)
     assert out_msg["data2"].shape == (3,)
     assert out_msg["data3"].shape == (3,4)

def _mgs(a):
    m, n = a.shape
    v = a.copy()
    r = numx.zeros((n, n))
    for i in range(n):
        r[i, i] = numx_linalg.norm(v[:, i])
        v[:, i] = v[:, i]/r[i, i]
        for j in range(i+1, n):
            r[i, j] = mult(v[:, i], v[:, j])
            v[:, j] = v[:, j] - r[i, j]*v[:, i]
    # q is v
    return v, r

 def test_inverse_message_routing(self):
     """Test the inverse routing for messages."""
     sboard = BiSwitchboard(input_dim=3, connections=[2,0,1])
     x = n.array([[1,2,3],[4,5,6]])
     msg = {
         "string": "blabla",
         "method": "inverse",
         "list": [1,2],
         "data": x,  # should be mapped by switchboard
         "data2": n.zeros(3),  # should not be modified
         "data3": n.zeros((3,4)),  # should not be modified
         "target": "test"
     }
     y, out_msg, target = sboard.execute(None, msg)
     assert y is None
     assert target == "test"
     reference_y = n.array([[2,3,1],[5,6,4]])
     assert out_msg["string"] == msg["string"]
     assert out_msg["list"] == msg["list"]
     assert (out_msg["data"] == reference_y).all()
     assert out_msg["data2"].shape == (3,)
     assert out_msg["data3"].shape == (3,4)

def _mgs(a):
    """Modified Gram-Schmidt."""
    m, n = a.shape
    v = a.copy()
    r = numx.zeros((n, n))
    for i in range(n):
        r[i, i] = numx_linalg.norm(v[:, i])
        v[:, i] = old_div(v[:, i],r[i, i])
        for j in range(i+1, n):
            r[i, j] = mult(v[:, i], v[:, j])
            v[:, j] = v[:, j] - r[i, j]*v[:, i]
    # q is v
    return v, r